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Preface 


A  common  problem  in  estimating  spares  requirements  for  complex 
systems  is  accurately  forecasting  failure  rates.  While  a  point  estimate 
of  system  reliability  can  normally  be  obtained,  it  is  relatively 
meaningless  without  some  confidence  limit  on  it.  The  problem  is 
compounded  in  highly  reliable  systems  that  have  some  period  of  failure 
free  life  before  entering  their  failure  period.  Considerable  work  has 
been  done  on  reliability  estimation  over  the  last  twenty  years  at  the  Air 
Force  Institute  of  Technology.  This  thesis  is  a  continuation  of  many 
efforts  to  provide  a  flexible,  robust  method  of  estimating  confidence 
limits.  In  an  attempt  to  make  the  methods  and  procedures  developed  as 
useful  as  possible  to  future  readers,  I  have  included  several 
illustrations,  an  example  of  how  to  use  the  procedures,  and  the  computer 
printouts  of  the  programs  used. 

Thanks  are  due  to  Dr.  Albert  H.  Moore  for  his  assistance  in 
selecting  this  topic  and  his  guidance  throughout.  Also,  the  frequent 
support  of  the  AFIT  and  ASD  computer  personnel  in  extending  my  account 
and  providing  timely  advice  on  methods  to  improve  turn-around  time  is 
much  appreciated. 


Murray  R.  MacDonald 
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A  Double  Monte  Carlo  method  of  obtaining  confidence  limits  for 
complex  systems  based  on  component  failure  data  assuming  a  three 
parameter  Weibull  distribution  was  developed.  Three  new  parameter 
estimation  routines  were  developed  and  compared  with  the  Harter-Moore 
three  parameter  maximum  likelihood  routine  for  use  with  the  Monte  Carlo 
method.  The  sensitivity  of  the  method  to  system  reliability,  sample 
size,  and  number  of  points  in  the  component  reliability  distributions  was 
assessed.  An  approximate  method  of  calculating  and  correcting  for 
parameter  estimation  bias  was  developed  and  illustrated.  The  Double 
Monte  Carlo  method  appears  to  be  effective  at  system  reliabilities  from 
74%  to  96%  with  componenent  failure  sample  sizes  as  small  as  five  with 
the  Linear  Least  Squares  parameter  estimation  routine  developed. 
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A  MONTH  CARLO  TECHNIQUE  SUITABLE  FOR  OBTAINING  COMPLEX  SPACE  SYSTEM 


RELIABILITY  CONFIDENCE  LIMITS  FROM  COMPONENT  TEST  DATA 
WITH  THREE  UNKNOWN  PARAMETERS 


I  Introduction 


Problem  Statement 

Accurate  estimation  of  complex  system  reliability  is  important  for 
operational  planning  and  system  replacement  scheduling.  Estimation  of 
system  reliability  for  Space  systems  is  particularly  important  in  that 
replacements  or  repairs  are  generally  difficult  to  make.  Little  system 
reliability  data  exists  due  to  the  small  numbers  of  each  type  of  system 
and  the  varying  complexity  of  the  systems.  Component  reliability 
continues  to  improve  and  longer  missions  arc  planned  and  conducted  (Ref 
32)  so  empirical  estimation  of  system  reliability  is  impractical,  thus 
necessitating  use  of  an  analytic  model  to  estimate  system  reliability. 
Most  models  simply  assume  a  constant  failure  rate  (exponential 
distribution)  for  the  components  (Ref  19),  but  these  models  have  proved 
to  be  very  conservative  which  results  in  larger  system  purchases  than 
required  (Ref  2).  Empirical  data  was  added  to  the  models  in  an  effort  to 
improve  their  accuracy  and  the  effect  of  mission  controller  selection  of 
alternate  system  modes  (work-arounds)  has  also  been  considered  (Ref  25) 
to  produce  reliability  estimates  which  are  not  grossly  in  error. 

Space  system  failures  have  tended  to  occur  in  one  of  two  separate 
periods:  early  system  failure  caused  by  undetected  defects,  and  wearout 
or  failures  caused  by  random  mishaps  later  in  the  system  life  (Ref  2:6). 
Accurately  predicting  early  failures  caused  by  defects  is  a  function  of 


quality  control.  With  improving  design  and  quality  control,  complex 
spacecraft  with  failure  modes  dominated  by  wearout  can  be  manufactured 
(Kef  2:7),  and  the  effect  of  a  few  defects  can  be  nullified  by 
work-arounds.  The  estimate  of  increasing  importance  is  the  effect  of 
long  terra  failures.  To  accurately  model  these  effects,  the  system 
component  failure  distributions  must  include  location  parameters. 

Incorporating  location  parameters  into  models  which  use  an 
exponential  distribution  assumption  would  account  for  the  necessary 
period  before  the  long  terra  failure  period  was  entered.  Unfortunately, 
the  exponential  distribution  is  not  robust  in  that  departures  from  the 
distribution  can  result  in  large  errors  (Ref  11  ).  The  primary  use  of  the 
exponential  distribution  in  system  reliability  estimates  would  appear  to 
be  for  a  system  which  is  composed  of  many  components  which  are  changed  on 
failure  (Ref  14:235*237).  This  is  not  the  case  for  space  systems,  so  an 
alternate  set  of  assumptions  which  can  more  accurately  model  the  failure 
distributions  encountered  must  be  used. 

A  point  estimate  of  system  reliability  is  relatively  easy  to  obtain 
but  is  of  itself  of  little  value  without  confidence  bounds.  For  example, 
a  system  reliability  estimate  of  0.99  with  a  90#  lower  bound  of  0.95  is 
considerably  different  than  a  system  reliability  estimate  of  0.99  with  a 
90#  lower  bound  of  0.50.  The  system  reliability  estimates  required  are  a 
point  estimate  and  a  lower  limit  at  some  pre-selected  confidence  level. 

The  problem  is  to  develop  a  robust  model  that  can  account  for  some 
period  of  guaranteed  life  before  entering  the  wearout  period  for 
reliability  confidence  interval  estimates  of  complex  systems.  The  model 
must  be  able  to  incorporate  many  different  types  of  components  with 
different  failure  modes  and  guaranteed  lives. 
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Review  of  Applicable  Literature 


Over  the  last  twenty  years,  a  considerable  amount  of  research  has 
been  done  into  estimating  component  parameters  from  failure  data  and 
system  reliability  from  component  data,  Orkand  in  I960  (Ref  24:18) 
suggested  a  Monte  Carlo  method  for  the  estimation  of  the  lower  confidence 
limit  on  the  reliability  of  a  complex  network  of  components.  He 
suggested  this  method  as  a  general  solution  procedure  and  provided  a  more 
detailed  solution  for  the  case  where  the  sample  data  for  each  component 
are  in  binomial  form.  Rosenblatt  in  1962  discussed  the  problem  of 
determining  confidence  limits  for  the  reliability  of  complex  systems. 

She  suggested  that  estimation  using  simulation  was  "...  the  simplest 
and  most  generally  applicable  procedure  for  estimating  R  through  F 

(xx . Xk)...M  (Ref  29;1iy)  and  used  a  binomial  theoretical  treatment  of 

the  problem. 

In  1963,  Quayle  (Ref  27)  summarized  the  applicable  reliability 
theory  and  provided  some  preliminary  work  on  parameter  estimation  with 
his  method  of  using  order  statistics  to  estimate  the  scale  parameter  of 
the  Weibull  probability  density  function.  The  same  year  Bernhoff  (Ref  3) 
showed  that  adding  component  confidence  limits  to  obtain  system 
confidence  limits  was  erroneous  and  that  no  single  system  parameter  was 
appropriate  when  the  components  have  different  distributions.  He 
determined  that  "The  analytical  solution  becomes  impractical  when  the 
system  reliability  estimator  is  the  function  of  two  or  more  dissimilar 
mathematical  forms  and  mathematical  simulation  must  be  used"  (Ref  3:3) • 
For  confidence  limits  he  generated  and  used  a  step  cumulative 
distribution  function. 

Levy  in  1964  (Ref  16)  used  a  Monte  Carlo  Technique  to  obtain  system 


reliability  confidence  limits  from  component  failure  test  data  assuming 
that  the  components1  failure  distributions  were  two-parameter  Veibull 
(location  parameter  =0).  He  established  »  step  cumulative  distribution 
function  to  obtain  system  lower  reliability  confidence  limits.  His  work 
was  later  consolidated  and  published  (Ref  17).  Moore  in  1965  (Ref  25) 
extended  the  concept  of  using  the  Monte  Carlo  technique  for  obtaining 
system  reliability  confidence  limits  from  component  data  for  cases  where 
the  mathematical  model  for  the  underlying  failure  distributions  is  known, 
component  test  data  exists  to  estimate  the  parameters,  and  the 
distribution  of  the  estimators  of  the  parameters  is  unknown.  The  basic 
method  consisted  of  obtaining  a  sample  distribution  of  reliabilities  from 
which  an  approximate  confidence  interval,  or  limit,  can  be  obtained  at 
any  level  of  confidence. 

In  1967  Hahn  and  Shapiro  in  their  text  (Ref  10:Chap  7)  discussed  the 
problem  of  estimating  confidence  intervals  for  complex  systems.  The 
methods  developed  were  the  use  of  the  Central  Limit  Theorem  for  series 
systems  with  a  large  number  of  components,  the  generation  of  system 
moments,  and  the  Monte  Carlo  method.  They  favor  the  generation  of  system 
moments  for  relatively  simple  systems  but  prefer  the  Monte  Carlo  method 
for  ”...  highly  complex  situations  for  which  the  method  of  generation  of 
system  moments  becomes  too  difficult.”  (Ref  8:246) 

In  1972  Lannon  (Ref  15)  used  the  Monte  Carlo  method  to  approximate 
system  reliability  confidence  limits  assuming  the  components  had  failures 
characterized  by  two-parameter  Weibull  distributions  (location 
parameterO) .  In  1975  Boardman  and  Kendall  (Ref  6)  developed  a  method  of 
parameter  estimation  for  a  binomial  mixture  of  two  single  parameter 
exponential  distributions  under  the  assumption  of  two  possible  causes  of 
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failure,  each  with  a  single  parameter  exponential  distribution.  Their 
method  may  have  some  application  to  the  problem  of  estimating  spacecraft 
reliability  but  suffers  from  the  assumption  of  location  parameters5^,  the 
lack  of  robustness  of  the  exponential  distribution,  and  the  simplicity  of 
the  model. 

In  197b  Gatliffe  (Ref  9)  extended  the  use  of  the  Log-gamma  procedure 
for  estimating  system  reliability  from  series  only  arrangements  to 
series-parallel  arrangements.  His  method  does  not  require  either 
assumptions  about  failure  distribution  of  any  component,  or  equal  sample 
sizes.  His  results  are  good  for  highly  reliable  systems  but  the  accuracy 
is  unstable  from  configuration  to  configuration.  His  system  also  can  be 
very  conservative  in  that  it  generates  artificial  failures  when  no 
failures  were  observed. 

Bilikam  and  Moore  provided  two  practical  illustrations  of  the  use  of 
the  Monte  Carlo  method  to  estimate  mission  reliability  in  1977  and  1978. 

In  the  first  (Ref  4)  they  used  time-grouped  mission  equipment  failure 
data  where  the  exact  failure  times  were  unknown  although  the  equipment 
was  known  to  have  failed  during  a  mission  of  known  length.  In  the  second 
(Ref  5)  they  used  known  failure  times  of  one  type  of  aircraft  engine 
component.  Also  in  1978,  Snead  (Ref  30)  developed  a  univariate  method  of 
using  the  asymptotically  normal  property  of  R(t)  with  a  Monte  Carlo 
technique  to  estimate  system  reliability  confidence  limits. 

In  1979  Putz  (Ref  26)  used  the  univariate  Monte  Carlo  method  to 
estimate  lower  confidence  limits  of  system  reliability  based  on  component 
test  data.  He  assessed  the  sensitivity  of  the  method  to  the  asymtotic 
normality  assumptions  and  estimated  the  minimum  sample  size  required  for 


this  method.  His  method  is  the  most  effective  when  component  and  system 


reliabilities  are  low  (less  than  0.9)  and  sample  sizes  of  fifty  or  more 
are  available.  When  the  component  reliabilties  are  high  and/or  the 
sample  size  is  low,  the  distribution  of  R(t)  is  no  longer  nearly  normal 
which  can  result  in  significant  errors# 

In  1979  Rice  (Ref  2ti)  assumed  the  number  of  component  failures  was 
binoraially  distributed*  Using  the  asyratotic  normality  property  of  the 
binomial  distribution  (n>20),  he  developed  a  Monte  Carlo  me  hod  of 
estimating  lower  confidence  limits  on  system  reliability  with  component 
failure  data  input.  In  the  cases  where  no  failures  were  observed  he  used 
the  Gatliffe  method  of  generating  artificial  failures#  Also  in  1979, 
Antoon  (Ref  1 )  used  Monte  Carlo  analysis  to  find  empirically  the  standard 
deviation  of  reliability  of  a  system  whose  underlying  component 
distributions  were  two-parameter  Weibull  (location  parameter=0) .  He 
developed,  by  curve  fitting,  an  equation  for  computing  the  standard 
deviation  in  terms  of  reliability  and  sample  size. 

In  I960,  Johnston  (Ref  13)  used  a  Modified  Double  Monte  Carlo 
procedure  to  estimate  system  reliability  from  component  data  where  the 
component  failure  distributions  were  characterized  by  the  two-parameter 
Weibull  (location  parameter-0 ) .  He  used  the  bias  tables  published  by 
Thoman,  Bain  and  Antle  (Ref  31 )  to  correct  his  estimates  of  reliability 
and  obtained  reasonable  results,  although  the  results  are  difficult  to 
evaluate  fully  since  different  system  configurations  with  different 
reliabilities  were  used.  Also  in  1930  Moore,  Harter,  and  Snead  (Ref  21 ) 
compared  three  Monte  Carlo  techniques  for  obtaining  system  reliability 
confidence  limits;  the  bivariate,  the  univariate,  and  the  Double  Monte 
Carlo.  Their  conclusion  was  that  the  bivariate  tended  to  be  conservative 


and  the  univariate  asymtotic  optimistic  with  the  Double  Monte  Carlo  in 


between.  Depuy  (Ref  7)  compared  the  accuracy  of  two  Monte  Carlo 
simulation  techniques  of  finding  lower  system  reliability  confidence 
limits;  the  bivariate  and  the  univariate.  She  found  the  bivariate 
method  the  most  accurate  if  the  true  system  reliability  is  below  0.95* 
and  the  univariate  most  accurate  if  the  true  system  reliability  is 
greater  than  about  0.9!?  and  the  component  data  sample  size  is  less  then 
twenty . 

Finally,  in  his  review  of  reliability  growth  (Ref  33),  Vonloh 
discusses  the  use  of  the  Monte  Carlo  method  in  system  reliability  growth 
prediction  models.  The  entire  subject  area  of  reliability  growth  is 
applicable  to  new  developing  technology  and  the  flexibility  of  th  Monte 
Carlo  method  in  general  makes  it  useful  for  estimation  of  the  reliability 
of  systems  whose  parameters  may  be  changing. 

Model  Selection 

Two  methods  of  system  reliability  determination  warrant  further 
discussion:  the  method  of  moments  and  the  Monte  Carlo  method.  The 
method  of  moments  can  be  the  most  economical  aproach.  It  also  allows  the 
analysis  of  the  relative  importance  of  each  component  variable  by 
examination  and  does  not  require  any  assumptions  about  the  underlying 
component  distributions.  On  the  other  hand,  the  accuracy  of  the  results 
is  not  always  consistent  and  cannot  be  readily  analyzed.  Also,  the 
generation  of  system  moments  soon  becomes  unworkable  with  increasing 
system  complexity  (Ref  10:246,247)* 

The  Monte  Carlo  method  requires  that  an  assumption  regarding  the 
component  failure  distribution  be  made.  Also,  it  does  not  allow  for 
detection  of  dominant  components.  Since  the  method  estimates  overall 
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system  performance,  a  change  in  any  of  the  components  requires  that  the 
entire  system  be  re-analyzed.  The  method  al30  requires  a  considerable 
amount  of  computer  time  -  the  exact  amount  being  dependent  on  the  system 
and  the  assumptions  in  the  model.  On  the  positive  side,  the  Monte  Carlo 
method  has  proven  useful  in  a  wide  variety  of  applications.  It  has  been 
extensively  used  in  developing  system  reliability  confidence  limits, 
particularly  for  a  two-parameter  Weibull  or  a  binomial  distribution.  It 
is  easy  to  use  and,  if  the  Double  Monte  Carlo  method  is  used  (Ref  20) 
does  not  require  any  assumptions  other  than  those  of  the  component 
reliability  distributions.  Therefore,  for  this  problem,  the  Double  Monte 
Carlo  method  was  selected  as  the  most  suitable  for  developing  the 
required  model. 

Underlying  Failure  Distribution  Selection 

In  using  the  Monte  Carlo  method,  it  is  necessary  to  select  a 
suitable  component  failure  distribution.  The  distribution  should  allow 
for  the  possibility  of  a  location  parameter  greater  than  zero  and  a 
non-symmetrical  shape  in  order  to  allow  fitting  of  the  distribution  to 
the  data  available  or,  as  Easterling  wrote  "Thus  the  task  facing  the 
statistician  is  more  often  one  of  model  fitting  than  of  parameter 
estimation"  (Ref  8).  This  is  because  a  given  set  of  data  may  not  clearly 
resolve  the  appropriate  distribution. 

Exponential.  The  exponential  distribution  is  widely  used  and  is 
well  analyzed.  However,  even  if  a  location  parameter  were  added  it  would 
still  not  have  the  required  flexibility  in  shape. 

Gamma.  The  gamma  distribution  has  been  used  in  fatigue  and  wearout 
studies.  It  can  assume  a  variety  of  forms  which  could  be  fitted  to  a 


considerable  variation  on  data. 


Normal.  The  normal  density  also  is  commonly  used.  It  can 
accomodate  a  period  of  near-guaranteed  life  and  can  assume  different 
scales  depending  on  the  mean  and  variance.  However,  its  symmetrical 
shape  limits  its  application. 

Log  Normal.  Like  the  normal  density  the  log  normal  can  accommodate 
a  period  of  near-guaranteed  life  and  can  assume  different  scales. 

However,  its  shape  is  limited  to  a  positively-skewed  normal  curve. 

Weibull .  The  three-parameter  Weibull  can  accomodate  any  positive 
location  parameter  and  a  wide  range  of  shapes  and  scale  depending  on  the 
respective  parameters.  The  scale  parameter  determines  the  spread  about 
the  mean,  the  shape  parameter  determines  the  failure  rate  -  whether 
increasing,  decreasing,  or  constant,  and  the  location  parameter 
determines  the  point  beyond  which  failure  can  occur.  The  Weibull  density 
function  has  shapes  that  are  similar  to  the  Gamma  or  the  lognormal  - 
assuming  appropriate  Weibull  parameters.  If  the  shape  parameter  is  1, 
the  Weibull  becomes  an  exponential  function;  a  shape  parameter  of  about 
3.7  yields  an  excellent  approximation  to  the  normal  function  and  a  shape 
parameter  of  2  can  approximate  Beta  distributions.  It  has  also  been 
shown  valid  for  a  wide  variety  of  actual  situations  (Ref  34)  and  has  the 
necessary  flexibility  to  fit  any  foreseeable  set  of  failure  data. 
Therefore,  the  three-parameter  Weibull  distribution  was  selected  for  this 
model • 

Objectives 

The  objectives  of  this  thesis  are: 

1  •  to  develop  a  model  to  estimate  complex  system  lower 


reliability  confidence  limits; 


2.  to  estimate  the  minimum  practical  sample  size,  and 

3.  to  assess  the  sensitivity  of  the  Double  Monte  Carlo  method 
to  the  number  of  points  in  the  sample  distributions  of  reliability. 

Assumptions 

It  is  assumed  that: 

1 .  The  underlying  component  failure  distributions  can  be 
modeled  by  three-parameter  Weibull  distributions; 

2.  Components  fail  independently;  that  is,  there  are  no 
secondary  failures; 

3*  A  mathematical  relationship  between  component  reliabilities 
and  system  reliabilities  can  be  established; 

4.  The  International  Mathematics  and  Statistics  Library  (IMSL) 
subroutines  GGUBS  and  GGWIB  provide  valid  random  variables;  and 

3.  The  user  has  a  basic  knowledge  of  reliability  theory,  Monte 
Carlo  methods,  and  FORTRAN  77- 

Approach 

Existing  methods  of  estimating  parameters  from  failure  data  were 
examined  and  the  most  suitable  method  was  selected  for  inclusion  in  the 
model.  The  failure  data  was  generated  artificially  to  represent  true 
component  failures  from  three-parameter  Weibull  distributions  with 
different  parameters.  A  single  complex  model  was  developed  and  the  true 
reliability  calculated  analytically  to  use  as  a  test  of  the  model 
results.  The  model  was  tested  at  system  reliabilities  of  about  75$,  05$, 
and  95$  with  all  of  the  component  reliabilities  roughly  matched  to 
simulate  a  balanced  system.  At  each  reliability  level  component  failure 
sample  sizes  ranging  from  five  to  fifty  were  modeled  to  assess  the  effect 
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of  different  sample  sizes  on  the  model  accuracy. 

The  Double  Monte  Carlo  method  was  used  to  generate  the  estimates  of 
system  reliability  confidence  intervals.  When  initial  results  indicated 
that  the  method  of  parameter  estimation  used  was  inadequate  for  a  three 
parameter  model,  three  new  methods  were  developed,  evaluated,  and  the 
best  selected  for  use.  The  evaluation  of  the  overall  method  consists  of 
comparing  the  percent  of  times  the  X  percent  confidence  interval  captures 
the  true  reliability.  For  example,  at  the  80^  confidence  level,  of 

the  time  the  confidence  level  should  be  below  the  true  reliability. 
Finally,  a  sample  illustration  of  the  method  was  provided  for  practical 
guidance. 
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Method  of  Maximum  Likelihood 

The  assumption  of  a  component  failure  distribution  necessary  for  the 
Monte  Carlo  method  requires  that  the  distribution  parameters  for  each 
component  be  estimated  from  its  failure  data.  The  method  of  maximum 
likelihood  has  been  widely  accepted  as  one  of  the  most  reliable  methods 
of  estimating  distribution  parameters.  The  maximum  likelihood  estimators 
are  consistent,  asymptotically  normal  and  asymptotically  efficient  for 
large  samples  under  most  conditions  (Ref  35:89). 

The  probability  density  function  of  a  random  variable  T  having  a 

Weibull  distribution  with  location  parameter  c,  scale  parameter  e,  and 

shape  parameter  k  is 

k  t  t-c  k 

f(t;c,e,k)  *  -  (“)  exP  [-  (■— )  ]  o<c<t,  e>o,  k>o 

\7  v  © 

To  establish  the  maximum  likelihood  values  of  the  parameters  c,  e, 
and  k  it  is  necessary  to  formulate  the  likelihood  function  and  solve  for 


the  values  of  the  parameters  that  maximize  the  function.  Let  ,  T2 > 
•••>  Tn  be  the  observed  values  in  a  random  sample  of  size  n.  Then  the 
likelihood  function  is 

L(t,c,e,k)  =  f(ti5  c,e,k)  tj_>c 

Now  if  the  t^  are  treated  as  fixed  constants,  then  the  likelihood 
function  may  be  treated  as  a  function  of  the  three  unknown  parameters. 
Substituting  in  the  value  of  f(t c,  e,  k)  the  likelihood  function 


The  natural  logarithm  of  the  likelihood  function,  In  L,  is  easier  to 
work  with  and  does  not  result  in  any  loss  of  generality  since  the  maximum 


of  In  L  and  the  maximum  of  L  will  occur  at  the  same  values  of  c,  ef  and 
k.  The  first  partial  derivatives  of  In  L  with  respect  to  each  of  the 
variables  ( the  three  unknown  parameters)  are  set  equal  to  zero  and  solved 
simultaneously  to  yield  the  maximum  likelihood  values  of  the  parameters. 
The  analytic  solution  of  this  system  of  partial  differential  equations  is 
intractable  and  so  requires  an  iterative  computer  routine.  A  commonly 
used  routine  which  has  proven  satisfactory  for  similar  applications  is 
the  Harter-Moore  method  of  false  position  (Ref  12).  This  routine 
estimates  the  maximum  likelihood  parameters  based  on  the  first  m  order 
statistics  of  a  sample  of  size  n  with  r  censored  from  below.  The 
formulation  of  the  natural  logarithm  of  the  likelihood  function  used  is 
ln(L)  *  ln(n!)  -  ln((n-m)!)  -  ln(r!)  ♦  (m-r) ( ln(k)-k  ln(e)) 

+  (k_l)  ilr,lln(ti_c)  'iir+1  Uti-c)/«]k 
-  (n-m)  L(tm-c)/eJk  ♦  rln  |l-exp  [-(tr*i -c)k/ek]} 


This  formulation  leads  to  the  partial  differential  equations 

3ln  L  k(m-r)  .  m  .  ,k  k  +  1  k(n-m)(tm-c) 

-  =  -  -  +  k.I  ,  (t.-c)  e  +  - : — 7^ - 

e  e  i~r+l  1  ek+l 

-  kr(t^+1-c)k  x  exp[-(  t^^1“-c)^/e^]/e^  +  1{  l-exp[-(tr+1-c)^/e^]} 

--"L  a  (m-r)  ( 1/k-lne)  +.?  ln(t.-c)-.?  t  [ (x.-cj/e]*1  ln[(x.-c)/e]  - 

3k  i  =r+ 1  i  i=r+l  i  l 

(n-m) [ ( t  -c)/e]k  ln[(t  -c)/e]  +  r(t  -c)k  ln[  (t  -c)/e] 
m  m  r+i  r+1 

exp{-f (tr+1~c) /e]k)+  ©k{ l-exp[-( tr+1~c)k/ek]} 


=  ( 1-k) .?  .  (t.-c)  1  +  k®  k  .?  (t.-c)k_1 

3c  i=r+l  i  i=r+l  l 

+  (n-m)ke  k  (t  -c)k  1  -  kr(t  -c)k  1  x 
m  r+i 

exp[-( trfl”C)k/ekJ  ek  {l-exp[-( tp+1-c )k/ek]  } 

The  routine  which  solves  these  equations  is  listed  in  Appendix  B. 
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Median  Rank  Values 


The  Double  Monte  Carlo  method  builds  estimated  reliability 
distribution  functions  from  which  the  appropriate  confidence  limits  are 
selected.  This  requires  that  the  randomly  generated  reliabilities  be 
ordered  and  ranked.  Several  methods  of  rank  plotting  are  available,  with 
the  median  rank  method  the  most  commonly  used  because  of  the  assumption 
that  the  rank  distributions  are  skewed.  The  median  rank  is  actually  an 
incomplete  beta  ratio  which  cannot  be  readily  calculated.  However,  the 
approximation  to  the  median  rank  value  given  by 


has  an  insignificant  error  for  the  large  sample  sizes  (n  >  50)  used  in 
the  reliability  distribution  functions.  An  illustration  of  the  median 
rank  plotting  against  reliability  with  linear  interpolation  between 
points  as  used  in  this  development  is  provided  in  Fig  1  . 


Figure  1 .  Median  Ranks 

Random  Variable  Generation 

The  IMSL  routine  GGWIB  was  used  to  generate  single  parameter  (k) 
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Weibull  variables  which  were  transformed  to  three-parameter  Weibull 
variables  by  subroutine  WhilBL.  The  IMSL  routine  uses  a  c.d.f. 

F( t)=1 -exp  t-tk]  and  inverts  this  to  the  reliability  function  to  obtain 
the  relationship  t=  [-In  (u)  )  where  u  is  a  uniform  pseudo-random 
variable. 

For  the  three-parameter  distribution,  the  reliability  function  is 

R(t)  =  exp[-  (-^)k] 

Letting  u=R(t)  and  taking  the  logarithm  of  both  sides  gives 

In  u  =  -(^)k 

(-In  u)^/k  =  — — 
e 

t  *  e(-ln  u)  1 +  c 


Double  Monte  Carlo  Method 

The  Double  Monte  Carlo  method  does  not  require  any  asymtotic 
assumptions  and  can  be  used  with  any  component  failure  data  providing 
that  an  assumption  is  made  regarding  the  underlying  failure  distribution. 
For  the  purpose  of  developing  and  proving  the  model,  the  "real"  failure 
data  was  generated  using  the  "true"  parameter  values  for  each  component. 
These  true  parameter  values  also  allowed  the  analytical  calculation  of 
the  true  reliability  which  was  used  as  a  test  of  the  results.  The  Double 
Monte  Carlo  method  initially  used  consisted  of  the  following  steps. 

1.  Generate  the  true  component  failure  data. 

2.  From  the  true  component  failure  data,  estimate  the  three 

parameters  of  each  of  the  component  reliability  functions. 

3.  Generate  a  simulated  sample  of  component  failures,  using  the 


15 


estimated  parameters,  with  the  same  number  of  observations  as  the 
test  data. 


4-  From  these  simulated  failures,  estimate  the  three  parameters 
of  each  of  the  component  reliability  functions. 

5*  Using  the  second  estimate  of  the  parameters,  calculate  the 
reliability  of  each  component, 

6.  Kepeat  steps  3-5  until  the  required  underlying  sample  size  is 
obtained  (50,  75  or  150). 

7*  Establish  sample  estimated  reliability  distribution  functions 
for  each  component  by  ordering  the  Rjj  for  each  component  and 
matching  each  with  the  appropriate  median  rank.  The  first  and 
last  order  statistic,  associated  with  the  median  ranks  0  and  1 
respectively,  are  approximated  using  linear  extrapolation  off  the 
two  nearest  order  statistics. 

8.  Randomly  select  a  reliability  for  each  component  from  its 
reliability  distribution  function  using  linear  interpolation  between 
points  and  compute  the  system  estimated  reliability  Repeat 

until  600  estimates  of  system  reliability  are  obtained. 

9*  Order  the  against  median  ranks  and  determine  the  99,  95,  90, 
80,  70,  60,  and  50  percent  lower  confidence  points  using  linear 
interpolation  between  points  on  the  system  sample  distribution  of 
reliability  estimates.  Note  if  the  true  reliability  is  greater  than 
or  equal  to  each  of  these  confidence  points. 

10.  Steps  1-9  provide  one  estimate  of  the  system  reliability 
confidence  limits.  To  validate  the  method,  these  steps  are  repeated 
1000  times.  The  X  percent  confidence  limit  should  be  less  than  or 
equal  to  the  real  system  reliability  X  percent  of  the  time. 
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Appendix  A  contains  the  computer  program  listing  for  the  Double 
Monte  Carlo  method  used.  In  actual  practice  only  one  true  set  of  data 
would  be  a  ailable;  step  9  would  consist  of  printing  out  the  desired 
confidence  levels  and  step  10  would  not  be  applicable. 


Ill  Method  Development 


System  Configuration  and  Reliability 

The  Monte  Carlo  method  of  determining  approximate  system  reliability 
confidence  limits  was  selected  as  the  most  suitable  for  this  development 
because  of  its  advantages  in  handling  complex  systems*  Of  interest  is 
the  accuracy  of  the  method  for  different  reliability  levels  and  different 
sample  sizes  for  a  complex  system.  Therefore,  only  one  system 
configuration  was  used  in  order  to  be  able  to  compare  the  results  at  the 
different  reliability  levels  and  sample  failure  levels  tested.  The 
system  configuration  selected  is  illustrated  in  Fig  2. 


Fig  2.  System  Configuration 

For  the  Weibull  distribution,  the  component  reliabilities  are 

f*i(t)  =  expt-C^^-1  )ki]  t>o  Ci>t  e^>o  k^>o 
ei 

From  the  system  illustration,  the  system  reliability  Rg  is 
Rs  =  R1  (1-FbFc)  =  R,  (1-(1-RB)(1-RC)  )  =■  R,  (RB  +  RC  -  RfiRc) 
fiB  *  RAR4 

Ra  =  (l-F^F-^)  ~  R2  +  “  R^R^ 


18 


RC  "  +  R6  -  V6 

Rs  =  B,  il (H^  +  U-}  -  H3)  R4]  +  (K5  ♦  R6  -  R5R6)  ~ 

L(R2  ♦  H3  -  «2  R5)  R4]  (H5  ♦  R6  -  R^f 

Parameters  and  Reliabilities  Selected 

As  a  matter  of  convenience  and  to  maintain  continuity  with  previous 
methods,  a  time  of  100  units  (l?  =  100)  was  used  throughout.  The  same 
location  parameters  were  used  with  the  scale  and  shape  parameters  changed 
to  provide  balanced  component  reliabilities  and  a  good  range  of  parameter 
selections  for  the  test.  The  parameters  and  reliabilities  used  are 
listed  in  Table  1. 

Estimates  Generated 

The  parameter  estimates  generated  by  the  method  of  maximum 
likelihood  are  biased,  and  the  reliability  estimates  derived  from  these 
estimated  parameters  will  be  biased .  Thozaan,  Bain,  and  Antle  (Ref  31) 
empirically  determined  and  tabled  the  bias  in  R(t)  for  a  two-parameter 
Weibull  distribution  for  a  range  of  .30  i  R(t)  1  .98  and  sample  sizes 
from  8  to  100.  Kor  sample  sizes  greater  than  15,  the  biases  were  only 
third  decimal  place  values.  Moore,  Harter,  and  Antoon  (Ref  22)  assess 
the  reliability  estimate  for  a  two-parameter  Weibull  as  being  very  nearly 
normal  and  very  nearly  unbiased  for  sample  sizes  greater  than  about  20. 
The  bias  in  R(t)  for  a  three-parameter  Weibull  distribution  has  not  been 
tabulated  but,  from  the  work  done  on  the  two-parameter  Weibull 
distribution,  can  be  expected  to  be  small  for  larger  sample  sizes  (>  30). 
The  results  of  this  method  include  any  bias  present,  and  should  provide 
some  feel  for  the  magnitude  of  the  bias  in  the  parameter  estimation 
routine  used. 
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TABLE  I 


Parameters  and  Reliabilities 


i  * 

1 

2 

3 

4 

5 

6 

Cj^  (Location) 

10 

0 

15 

30 

25 

n 

(Shape) 

2.8 

1  .1 

2.0 

2.5 

1  .7 

mSM 

(Scale) 

140 

400 

180 

120 

180 

1{S  «  0.741 

.748 

.804 

.800 

.771 

.798 

ci 

10 

0 

15 

30 

25 

50 

ki 

2.5 

0.8 

2.0 

3.5 

1  .2 

2.0 

•i 

200 

470 

180 

120 

250 

150 

Ri 

R3  =  0.849 

.853 

.748 

.800 

.859 

.790 

.895 

ci 

10 

0 

15 

1 

25 

50 

ki 

2.9 

2.1 

2.3 

■ 

2.7 

2.0 

•i 

270 

470 

400 

■ 

250 

280  j 

Rs  =  0.959 

.960 

.962 

.972 

■ 

.962 

.969 

1 

Accuracy  of  Method 

With  Monte  Carlo  simulation,  the  larger  the  number  of  trials  in  the 
simulation,  the  more  precise  the  solution  will  be*  The  desired  degree  of 
precision  can  be  obtained  by  increasing  the  number  of  trials.  The  number 
of  trials  required  for  a  degree  of  precision,  E,  can  be  calculated  at  a 
desired  confidence  level,  1  ,  by  considering  the  Monte  Carlo  as  a 

binomial  problem  where  the  estimate  of  interest  is  the  proportion  p  of 
systems  above  a  certain  level.  The  calculation  of  n  depends  on  the  value 
of  p  actually  found  in  the  Monte  Carlo  simulation  so  a  certain  amount  of 
trial  and  error  is  required.  However,  for  a  conservative  estimate  the 
largest  n  will  be  required  for  p  =  0.5  which  may  be  used  to  obtain  an 


upper  bound  on  n.  The  normal  approximation  to  the  binomial  can  be  used 
for  convenience  with  little  error  since  the  number  of  Monte  Carlo  trials 
will  generally  be  sufficiently  large  to  ensure  that  np  or  n  (l-p)  are 


greater  than  five.  This  approximation  leads  to  the  relation 


n 


PU-P*  z2 

1 -«/2 


for  a  two  sided  interval  where  is  the  (l-<»/2)l00  percent  point  of 

the  standard  normal  distribution* 

If  the  error  is  required,  it  can  be  determined  by  the  relation 


K 


z 


1  -«/2 


For  example,  if  1000  Monte  Carlo  trials  result  in  900  points  within  some 
specified  tolerance,  p  =  *  0*9  and,  at  the  35%  confidence  level 


E 


(.9)(.l)  % 
1000 


1.96 


±  0.0186 


While  the  accuracy  of  the  Monte  Carlo  procedure  can  be  readily 


estimated,  the  model  accuracy  is  more  dependent  on  the  accuracy  with 
which  the  component  parameters  can  be  estimated  from  failure  data*  An 
error  in  parameter  estimation  is  compounded  by  the  use  of  these  estimated 
parameters  to  generate  random  samples  from  which  the  second  estimate  of 
parameters  is  made.  Therefore,  the  overall  accuracy  of  the  method  can 


only  be  estimated  from  the  results  tested  against  the  known  point.  An 


underlying  assumption  in  developing  the  model  is  that  the  results  can  be 
extended  to  other  similar  complex  systems. 
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IV  Preliminary  Results  and  Method  Development 
Preliminary  Results 

The  Double  Monte  Carlo  program  was  developed  and  run  for  a  single 

estimate  of  system  reliability  (in  step  10  only  one  estimate  obtained)  at 

*  .95  with  10  failures  and  75  points  in  the  sample  distributions  of 

component  reliability  estimates.  This  single  estimate  required  677 

seconds  of  CDC  L^yber  Model  74  (CSB  System)  CPU  time.  When  extrapolated 
* 

to  1000  estimates  ‘per  run  and  15  runs  required  (5  sample  sizes  and  3 
reliabilities),  this  results  in  an  estimate  of  over  2820  hours  (4  months) 
of  CPU  time.  The  expensive  part  of  the  method  was  step  Four: 
calculating  the  second  maximum  likelihood  estimators  of  the  parameters 
from  the  simulated  failures.  In  addition  to  being  time  consuming,  this 
step  produced  estimates  of  the  location  parameter,  c,  that  were  larger 
than  the  test  time  (T  *  100)  for  202  of  the  450  (45$)  parameter 
estimations  from  the  simulated  failure  data.  These  large  estimates  of  c 
were  not  surprising  in  light  of  the  high  component  reliabilities  and  the 
small  sample  size  but  did  indicate  a  potential  problem  with  the  parameter 
estimation  method  under  these  circumstances. 

Parameter  Estimation  Development 

These  initial  results  showed  that  a  much  faster  and  more  reliable 
method  of  parameter  estimation  had  to  be  obtained  in  order  to  reduce  the 
CPU  time  and  get  more  reliable  estimates  of  the  location  parameter  c. 
Because  of  the  large  number  of  parameter  estimations  required,  the 
overriding  requirement  was  to  greatly  increase  speed.  Other  work  done  on 
parameter  estimation  was  more  concerned  with  accuracy  (Ref  18),  so  the 
other  available  routines  were  also  slow.  The  first  approach  taken  was  to 
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modify  the  maximum  likelihood  routine  used.  Initially  the  accuracy 
tolerances  were  set  at  .0001  and  the  program  could  run  for  a  maximum  of 
550  iterations.  Trial  and  error  with  the  method  resulted  in  reducing  the 
maximum  number  of  iterations  to  300  and  the  accuracy  tolerances  to  .01 
without  any  significant  degradation  of  results.  This  reduced  the  run 
time  to  about  60$  of  its  previous  level.  Since  this  was  still  far  too 
slow,  modifying  this  procedure  was  abandoned  and  three  new  methods  were 
developed . 

The  first  approach  taken  was  to  develop  a  computerized  graphical 
estimation  technique  using  the  ordered  samples  t^»  i  =  1f  2,  . ..,  n,  and 
accept  the  parameters  which  gave  the  minimum  error  least  squares  fit. 

The  cumulative  distribution  function  for  the  two-parameter  Weibull, 
F(t)  *  1-exp[-  (-)K]  can  be  rearranged  and  the  natural  logarithm  taken 
twice  to  give  the  relationship 

In  ( ln( fZp ('t]  ) )  =  k  ln(t)  "  k  ln(«) 
which,  when  rearranged,  gives 

ln(t)  -  jj  ln(ln(i:--~))  ♦  ln(a) 

The  substitution  of  Y  *  ln(t),  m=^ ,  X  =  ln  ( ln(  ) ) ,  and  a  *  ln(e) 

provides  the  linear  relationship  Y  *  mX  ♦  a.  The  values  of  F(t)  were 
estimated  by  the  use  of  median  ranks.  The  value  of  k  was  estimated  from 
the  first  and  last  values  on  the  abscissa  and  ordinate.  Using  the 
estimate  of  k  and  F(t)  as  given  constants,  ln(e)  can  be  calculated  by 
ln(«)  =  ln(t)  -  £  In  (ln(— ~y)) 

Since  e  is  a  constant  for  each  set  of  data,  the  value  of  ln(e)  is 
constant  for  each  t^  which  leads  to  the  following: 


it,  !«(•)  -J,  LMtj)  -  l  in 
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a  ln(e)  =  ln(ti)  -  £  In  (ln( 

ln(e)  =  (i?1  ln^)  -  ~  In  ( ln(  1_F,|“t“y ) )  >/n 

The  value  of  ln(e)  can  be  used  to  estimate  In  (t^)  and  these 
estimated  values  compared  to  the  actual  values  of  ln(t^)  observed  for  one 
estimated  linear  relationship.  If  various  values  of  c  are  subtracted 
from  the  sample  data,  each  new  sample  can  be  used  to  estimate  k  and  e  and 
the  parameters  which  provide  the  best  fit  to  a  linear  relationship 
accepted.  The  method  allows  for  an  estimate  of  c=t^  by  using  t2 
instead  of  tj  and  n  respectively  wherever  required.  In  this  case,  t^  is 
effectively  censored. 

The  step-by-step  procedure  for  this  method  is  as  follows: 

1.  Generate  the  abscissa  values,  x^»  from  In  ( 1 n ( ^ tf  )  ^  ^  u8*-n& 
median  ranks  as  the  plotting  position  of  F(t^). 


2. 

Calculate  the  ^Z^  x^ 

3. 

In  a 

loop  ranging  from  J  =  0  to  10  do  the  following: 

a. 

set  c j  -  ( J )  (0,1 )  ( t^) 

b. 

calculate  the  yi  by  the  relationship  y*  *  ln(t^-Cj) 

c . 

n 

calculate  • Z ^  y^ 

d. 

set  k*  =  2nl^l 

j  yn-yi 

e. 

calculate  In  (e)  by  ln(ej)  =  (i£1  yt  -  £  x±)/a 

f. 

calculate  the  natural  logarithm  of  the  estimated  failure 

times, 

B(yi), 

by  E(yi)  =  l-  *  ln(«j) 

«• 

calculate  the  sum  of  squares  of  errors,  SSEj,  by 

SSEj  •  ill  (E(yi)  _yi)2 

4.  Determine  the  smallest  SSEj  and  use  c*Cj,  k*kj,  e*exp( ln(ej ) ) 
The  method  of  estimating  k  from  the  slope  calculated  with  the 
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1 


extreme  values  could  result  in  some  degradation  of  the  estimate  if  one 
(or  both)  of  the  extreme  values  were  outliers.  Due  to  their  long 
separation,  a  considerable  deviation  would  have  to  occur  in  one  of  the 
extreme  values  of  t  before  the  effect  could  be  expected  to  be  serious. 
However,  to  avoid  any  potential  problem  of  this  nature,  or  others  pointed 
out  further  on,  it  is  only  prudent  to  examine  the  data  before  entering  it 
into  the  computer  for  analysis.  Because  of  the  method  of  obtaining 
linear  relationships  and  checking  them  by  best  least  squares  fit,  this 
method  will  be  referred  to  as  the  Linear  Least  Squares  (L.L.S.)  method  of 
parameter  estimation.  Appendix  B  contains  the  program  for  this  method  of 
parameter  estimation. 

Since  the  L.L.S.  method  depends  to  some  extent  on  the  plotting 
position  used,  various  plotting  positions  were  tested  against  several 
parameter  combinations  similar  to  those  selected  for  the  Double  Monte 
Carlo  procedure  developed.  The  plotting  positions  tested  were: 

1 .  the  mean,  j/ (n+1 ) ; 

2.  median,  ( j- .3)/(n+ .4 ); 

3.  (j-.375)/(n*.25); 

4*  midpoint,  (j-.5)/n;  and 
the  mode,  ( j-1 )/( n-1 ) . 

The  median  rank  method  was  selected  since  it  gave  consistently 
closer  estimates  than  the  other  two  methods  and  also  it  has  been 
extensively  used  with  good  results. 

Initial  LLS  Parameter  Estimation  Tests 

Three  manual  plots  of  twelve  data  points  each  on  Weibull  paper  were 
made  as  a  check  for  gross  errors  in  the  LLS  method.  For  the  first  check, 
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parameters  of  0=10,  k=2.3  and  e-120  were  used  to  generate  the  data 
points.  The  LLS  routine  provided  estimates  of  c=40.3>  k=1.48,  and  ©=105* 
The  plot  of  the  data  points  with  c=0  (•)  and  c=40.3  (©)  is  provided  in 
Figure  3.  The  estimate  selected  by  the  program  appears  to  plot 
accurately  on  the  Weibull  chart.  For  the  second  and  third  checks,  the 
data  points  corrected  for  c  which  the  program  selected  were  plotted  and 
the  parameters  estimated  manually.  In  both  of  these  cases,  the  data 
plotted  well  and  there  was  no  appreciable  difference  in  the  plotted 
parameter  estimates.  Parameters  used  to  generate  the  data  were  0*30, 
k=2.3»  ©=180,  and  c=30,  k=3.1  and  ©=400.  The  respective  parameter 
estimates  were  c=81 ,  k=1.35,  e=158  and  c=162,  k=1.76,  ©=302.  The  plots 
of  the  estimates  are  provided  in  Figures  4  and  5. 

The  LLS  method  of  parameter  estimation  was  then  tested  by  the 
following  computerized  procedure. 

1 .  Generate  random  samples  of  size  5  to  50  from  three-parameter 
Weibull  distributions  with  pre-selected  parameters. 

2 •  Calculate  true  system  reliability  from  the  pre-selected 
parameters.  T  was  arbitrarily  set  equal  to  100  throughout. 

3-  Estimate  the  parameters  from  the  sample  data  generated. 

4.  Calculate  the  estimated  reliability,  R(100),  from  the  estimated 
parameters.  If  c>100,  set  R(100)  =  1. 

5.  Three  sets  of  parameters  were  each  used  to  generate  5,  10,  20, 
and  50  random  samples  for  fifty  parameter  estimations  at  each  of  the 
twelve  combinations. 

6.  The  mean  square  error  between  the  estimated  reliabilities  and 
the  true  reliability  was  calculated  for  each  combination. 

The  true  parameters  used  for  the  test,  true  reliabilities,  number  of 
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time  R(100)=1,  and  mean  square  errors  for  reliability  are  listed  in  Table 
H.  A  graphical  illustration  of  the  method  from  ten  random  samples  is 
provided  in  Figure  6.  For  this  illustration,  the  samples  were  105.3, 
133.6,  196.4,  233.6,  360,  365,  372.1,  417.2,  428.5,  and  566.8  from  true 
parameters  of  c=25,  k=2.1,  and  e=280  with  a  true  reliability  of  0.939- 

The  parameters  obtained  from  the  '‘inear  least  squares  method  fit  the 
estimated  cumulative  distribution  function  to  the  true  cumulative 
distribution  function  very  well.  However,  the  estimated  parameter  sets 
generally  had  a  c  that  was  low,  and  a  k  and  e  that  were  high.  For 
example,  with  c=50,  k=2  and  e=280  and  a  sample  size  of  ten,  the  average 
estimated  parameters  were  c=3-84,  k=2.86  and  e=331 •  Figure  7  illustrates 
a  cumulative  distribution  obtained  from  the  true  parameters  and  the 
average  estimated  parameters.  This  tendency  to  underestimate  c  and 
overestimate  k  and  e  was  evident  for  all  the  initial  sample  size  and 
parameters  tested. 


Fig  6.  Illustration  of  Graphical  Method 
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Sample  Size 

MSE 

R ( 1  00 )  =  1 

c 

k 

© 

R(100) 

5 

.0234 

0 

10 

2.8 

140 

.7481 

10 

-0093 

0 

20 

.0063 

0 

50 

.0035 

0 

5 

.0206 

1 

0 

1  .1 

400 

.8044 

10 

.0083 

0 

20 

.0052 

0 

50 

.0032 

0 

5 

.0033 

3 

50 

2.0 

280 

.9686 

10 

.0009 

1 

20 

.0008 

22 

50 

.0002 

1 

0 

_ 1 
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As  a  further  test  of  the  L.L.S.  method,  a  comparison  of  individual 
parameter  MSEs  was  made  with  the  results  of  Miller  (Ref  18)  using  the 
parameters  c=10,  e=1  and  k  ranging  from  .5  to  4.0.  The  results  with 
these  parameters  were  poor  so  the  test  was  discontinued  for  further 
analysis  of  the  preliminary  findings.  An  attempt  was  made  to  plot  the 
data  points  on  Weibull  paper  but  this  proved  infeasible.  With  the  value 
of  c  extremely  high  compared  to  e,  all  of  the  points  effectively  plotted 
together.  Any  least  squares  errors  from  the  linear  relationship  under 
these  conditions  becomes  rather  meaningless.  When  the  data  is  run  on  a 
computer  program,  the  true  best  estimate  was  close  to  the  first  order 
statistic,  but  because  of  the  clustering  of  data  points  with  c  close  to 
0,  the  computer  would  often  select  a  low  value  of  c  with  a  resulting  high 
value  of  e.  Table  III  is  an  example  of  the  data  points  evaluated  on  one 
run  of  the  L.L.S.  program.  The  points  are  values  of  failure  times  (t) 
minus  c  over  a  range  in  c  from  c=0  to  c=t(l)  in  increments  of  1/10  of  the 
first  order  statistic,  t(l).  For  this  example,  the  following  values 
applied:  sample  size  =  12 

c  =  10 
k  »  2 
e  *  1 

MSE  c  =  62.9 
MSE  k  =  207 
MSE  e  -  66 
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Example 

of  Data 

for  L.L.S. 

.  Teat 

c=0 

10.34 

10.46 

10.56 

10.59 

10.67 

10.86 

10.93 

11.27 

11.29 

11.44 

11.5  2 

12.56 

c-.l  t(1  ) 

9-31 

9.43 

9.52 

9.56 

9-63 

9.82 

9-95 

10.23 

10.26 

10.41 

10.48 

11.53 

c=  .2  t(l  ) 

8.27 

8.39 

8.49 

8.52 

8.60 

8.79 

8.91 

9.20 

9.22 

9.38 

9.45 

10.49 

c=.3  t(l  ) 

7.24 

7.36 

7.45 

7.49 

7.56 

7.75 

7.88 

8.16 

8.19 

8.34 

8.41 

9.46 

O 

II 

. 

-P* 

c+ 

6.20 

6.32 

6.42 

6.45 

6.53 

6.72 

6.84 

7.13 

7.15 

7.31 

7.38 

8.43 

c=.5  t(l ) 

5.17 

5.29 

5-38 

5.42 

5.49 

5.68 

5.81 

6.10 

6.12 

6.27 

6.34 

7.39 

c= .6  t(  1  ) 

4.13 

4.25 

4.35 

4.38 

4.46 

4.65 

4.77 

5.06 

5.09 

5.24 

5.31 

6.36 

c=.7  t(l  ) 

3.10 

3 .22 

3.32 

3.35 

3.43 

3.62 

3.74 

4.03 

4.05 

4.20 

4.27 

5.32 

c=  .8  t(l) 

2.06 

2.18 

2.28 

2.31 

2.39 

2.58 

2.70 

2.99 

3.02 

3.17 

3-24 

4.29 

o 

II 

v£> 

oP 

1 .03 

1.15 

1 .25 

1 .28 

1 .36 

1 .55 

1  .67 

1 .96 

1 .98 

2.13 

2.21 

3.25 

c=t(l  ) 

.120 

.216 

.250 

.326 

.517 

.640 

.928 

.952 

1.10 

1.17 

2.22 

The  second  method  developed  was  actually  derived  from  the  L.L.S. 
method.  In  order  to  avoid  the  problems  caused  by  a  large  c  relative  to 
e,  c  was  estimated  from  a  linear  extrapolation  of  the  first  two  order 
statistics.  This  eliminated  the  time  consuming  process  of  selecting  the 
best  fit  to  a  linear  relationship,  so  k  was  estimated  by  the  average 
slope  between  order  statistics,  e  was  estimated  as  before  using  the 
values  of  c,  k,  and  the  order  statistics.  This  modified  L.L.S.  method 


proved  to  be  extremely  fast  in  the  parameter  estimation  comparison  tests 
described  later.  The  computer  routine  for  this  modified  L.L.S. 

(M.L.L.S.)  is  included  in  Appendix  B. 

The  third  new  method  of  parameter  estimation  developed  used  the  same 
procedure  of  extrapolating  from  the  first  two  order  statistics  to 
estimate  c.  Once  c  is  known,  k  and  e  can  be  quickly  estimated  by  maximum 
likelihood.  These  estimators  of  k  and  e  are: 

k  =  n/j  ([n^  Xj>  lnU^J/jSj  ln(xi)} 

•  =  Ujj  xik)/nJ1/k 

An  iterative  routine  was  used  to  find  k  where  x^^t^-c.  The  estimate 
of  k  was  then  used  directly  to  estimate  e.  This  method  is  referred  to  as 
the  Modified  Maximum  Likelihood  (M.M.L.)  method.  The  computer  routine 
for  the  M.M.L.  method  is  included  in  Appendix  B.  Theoretical  development 
of  these  estimators  of  k  and  e  is  included  in  Appendix  C. 

Parameter  Estimation  Tests 

The  modified  L.L.S.  method  and  the  M.M.L.  were  tested  for  individual 
parameter  MSEs  using  ca10,  e=1 ,  and  k  from  0.5  to  4.0.  The  results  of 
these  tests,  listed  in  Table  IV,  are  comparable  to  those  obtained  by 
Miller  (Ref  18)  for  these  parameter  selections.  However,  for  this 
project  the  accuracies  of  estimation  of  the  parameters  of  more  widely 
spread  distributions  were  required.  Therefore,  the  component  parameters 
previously  selected  for  an  855t  system  reliability  were  used  to  test  the 
L.L.S.,  the  modified  L.L.S.,  the  L.L.S.  using  an  average  slope  for  k 
instead  of  just  the  extreme  values,  and  the  M.M.L.  Sample  sizes  of  5, 

10,  20,  and  50  were  selected  to  provide  a  wide  range  without  an  excessive 
number  of  points.  All  of  the  above  tests  consisted  of  calculating  the 
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Sample  Size 


5  10  20  50 


Jj  .  L  .  S  . 

c 

959* 

906 

|  2152 

900* 

k 

11.4 

4.30 

3.83 

2.47 

0 

1175* 

1037* 

1 701  * 

944* 

L.L.S. 

c 

996 

905* 

2120 

900 

Average 

k 

9- S3 

3.15 

2.92 

1  .36 

Slope 

0 

1191 

1068 

1850 

1013 

2545 

1880 

9§er 

E 

3-47 

2.81 

1 .67 

0 

1  4455 

3184 

2375 

1231 

M.M.L. 

c 

2572 

mgBsm 

1016 

k 

2.80* 

1  .32* 

0 

3191 

1209 

c  »  50  k  *  3.5  9-120 


Sample  Size 


5  10  20  50 


L.L.S. 

c 

2428 

2402 

1385 

2500 

k 

5.16 

2.33 

.923 

2.25 

0 

3816*  I 

3350 

1125 

2686 

L  •  L .  S  • 

c 

2425* 

2322 

1264 

2500 

Average 

k 

4.43 

1  .63 

.539 

.959 

Slope 

0 

3931 

3441 

1083* 

3343 

c 

2998 

1545 

879 

294* 

L.L.S. 

k 

.763* 

.618 

.409 

.189 

0  ! 

4387 

2615 

1583 

570 

M.M.L. 

c 

2298 

i"536* 

813* 

309 

k 

.801 

.586* 

.319* 

.134* 

0 

4514 

2571* 

1429 

542* 

c  =  50  k  -  2.0  9-150 
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MSEs  based  on  Monte  Carlo  runs  of  1000  repetitions.  The  results  of  the 
estimations  are  listed  in  Table  V  with  the  best  parameter  estimate  for 
each  sample  size  marked  by  an  *.  The  computer  CPU  times  required  for 
1000  runs  of  the  four  sample  sizes  for  each  method  are  listed  in  Table 
VI. 


TABLE  VI 

Parameter  Test  CPU  Times 
Method  Time  in  Seconds 


L.L.S. 

76 

Modified  L.L.S. 

22 

L.L.S.  Average  Slope 

90 

M.M.L. 

46 

For  comparison,  a  limited  test  of  the  Harter-Moore  method  of  maximum 
likelihood  program  was  run.  On  the  initial  runs  using  1000  Monte  Carlo 
repetitions,  c=0,  k=0.8,  and  e=470,  only  the  runs  for  sample  sizes  5  and 
20  converged  within  400  seconds  so  no  further  runs  of  1000  repetitions 
were  made.  The  MSEs  from  these  two  runs  are  listed  in  Table  VII. 


TABLE  VII 

Harter-Moore  Method  -  1000  Repetitions 

Sample  Size  MSE  c  MSE  k  MSE  e 

5  9013  .3978  73574 

20  342  .025  17178 

The  number  of  repetition  was  decreased  to  100  and  a  further  test  was 
run  using  sample  sizes  of  5,  10,  20,  and  50,  c*25,  k=1.2,  and  e*250.  The 
run  at  sample  size  10  took  70  seconds  which  was  proportional  to  the 
length  of  time  taken  for  the  other  runs.  The  MSEs  of  this  test  are  listed 
in  Table  VIII. 
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TABLE  VIII 


Harter-Moore  Method  -  100  Repetitions 
Sample  Size 


3 

10 

20 

50 

MSE  c 

4100 

1238 

353 

156 

MSE  k 

1 .173 

.329 

.092 

.033 

MSE  o 

1  3333 

9947 

2735 

959 

c  =  25  k  =  1  .2  e  =  25 0 

Discussion  of  Parameter  Estimation  Test 

On  the  basis  of  the  tests  using  c=10  and  e=1  there  is  little 
difference  between  the  modified  L.L.S.  and  the  M.M.L.  in  accuracy,  with 
the  modified  L.L.S.  being  twice  as  fast.  The  estimates  of  c  are  the  same 
for  both,  with  the  occasional  small  differences  attributable  to  different 
random  number  seed  values,  since  they  use  the  same  method  of  estimating 
c.  Both  methods  give  reasonably  good  results  and  both  are  easy  to 
implement.  The  L.L.S.  method  does  not  produce  satisfactory  results  with 
a  large  c  relative  to  9.  if  the  data  is  all  clustered  about  one  point, 
then  some  method  must  be  used  to  spread  it  for  analysis  such  as 
estimating  a  location  parameter  from  the  first  two  order  statistics  or 
subtracting  a  large  fraction  of  the  first  order  statistic  from  all  the 
data. 

Using  data  generated  from  more  widely  spread  parameters,  the  L.L.S. 
appears  very  good.  It  is  more  consistent  in  its  estimates  of  c  and  the 
estimate  of  k  and  e  are  as  good  or  close  to  as  good  as  the  estimates  of 
any  of  the  other  methods.  It  even  compares  favorably  with  the  well 
established  Harter-Moore  method  of  three  parameter  maximum  likelihood 
estimation,  particularly  for  small  sample  sizes.  Of  interest  is  that  the 
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L.L.S.  generally  provided  better  estimates  than  the  L.L.S.  modified  to 
use  an  average  slope  to  estimate  k.  The  M.M.L.  method  was  more 
consistent  in  its  estimates  of  k  than  any  of  the  other  methods.  A  simple 
count  of  the  number  of  parameters  most  accurately  estimated  shows  the 
L.L.S.  and  the  M.M.L.  method  tied  at  25  each  with  the  L.L.S.  better  at 
estimating  c  and  the  M.M.L.  better  at  estimating  k.  Neither  parameter 
can  be  considered  the  more  critical  because  the  importance  of  an  error  in 
either  depends  on  their  magnitude  and  the  relative  magnitudes  of  c  and  e. 
The  selection  of  the  best  method  depends  on  which  produces  the  best 
results  in  the  Double  Monte  Carlo  procedure.  Given  equivalent  results, 
the  fastest  method,  M.M.L.,  would  be  preferable. 

The  results  of  the  tests  on  parameter  estimation  methods  indicate 
strongly  that  extensive  trials  with  a  wide  range  of  combinations  of  the 
three  parameters  are  necessary  before  any  method  can  be  selected  as  the 
universal  “best"  or  even  as  the  best  for  a  particular  purpose.  Until 
this  is  done,  it  is  advisable  to  visually  inspect  the  data  prior  to  any 
furhter  analysis.  If  the  data  appears  well  spread,  the  L.L.S.  would  be  a 
good  choice  of  method  to  estimate  the  parameters.  If  the  data  appears  to 
be  bunched  in  one  small  range,  then  the  M.M.L.  would  likely  be 
preferable.  Also,  as  for  any  other  use  of  data,  the  data  should  be 
checked  for  outliers  before  parameter  estimation,  and  for  goodness  of  fit 
after  parameter  estimation. 

The  original  systems  reliability  confidence  level  tests  were 
selected  for  sample  sizes  ranging  from  ten  to  one  hundred.  These  sizes 
were  selected  primarily  because  of  the  need  for  relatively  large  sample 
sizes  to  obtain  consistent  results  with  the  maximum  likelihood  method  of 
parameter  estimation.  However,  the  results  obtained  from  the  L.L.S.  and 
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the  M.M.L.  methods  developed  appear  reasonable  for  all  the  sample  sizes 
tested:  5,  10,  20  and  50.  Therefore,  the  project  was  modified  to  use 
these  sample  sizes  rather  than  the  sample  sizes  originally  picked. 

To  select  the  better  method  and  to  illustrate  the  use  of  the 
program,  1000  Monte  Carlo  runs  were  divided  into  seven  sets  of  "real" 
data  for  testing  at  each  of  the  system  reliabilities  74$,  85$,  and  96$. 
Each  of  these  sets  used  the  real  data  once  but  generated  simulated 
failures  and  confidence  limits  1000  times.  Also  checked  on  these  runs 
was  the  number  of  times  the  estimate  of  c  was  greater  than  the  selected 
mission  time  of  100  units.  For  the  L.L.S.,  the  largest  percent  of  times 
was  about  5$  compared  to  about  22$  for  the  M.M.L.  The  overall  results 
for  the  L.L.S.  method  were  much  better  than  the  results  for  the  M.M.L. 
method  at  85$  system  reliability,  so  the  M.M.L.  method  was  not  tested 
further.  For  this  project,  the  L.L.S.  method  of  parameter  estimation  was 
selected  as  the  best.  The  average  results  of  these  test  runs  are  listed 
in  Table  IX. 
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TABLE  I* 


Parameter  Estimation  Selection  Results 
L.L.S.  Method 


R  =  74# 

Confidence  Level  Sample  Size 


.50 

« 60 

.70 

.80 

.90 

.95 

.99 

5 

10 

20 

50 

.297 
.418 
.360 
•  345 

.387 

.499 

.396 

.478 

.662 

.460 

.443 

.653 

.853 

.622 

.640 

.782 

.975 

.793 

.885 

.851 

1 .0 

.899 

.983 

.911 

1 .0 

.993 

1 .0 
.995 

R  * 

85# 

Confidence 

Level 

Sample  Size 

.50 

.60 

.  70 

.80 

.90 

.95 

.99 

.313 

.428 

.729 

.904 

.998 

1  .0 

1  .0 

5 

•  445 

.473 

.628 

.850 

.962 

1  .0 

10 

.280 

.380 

.441 

.541 

.589 

.785 

.978 

20 

.335 

.598 

.707 

.804 

.875 

.975 

1  .0 

50 

R  * 

96# 

Confidence 

Level 

Sample  Size 

.50 

.60 

.70 

.80 

•  90 

.95 

.99 

•  345 

.490 

.759 

.977 

1 .0 

1 .0 

1 .0 

5 

.449 

.450 

.456 

.608 

.948 

.998 

1 .0 

10 

.197 

.249 

.251 

.268 

.511 

.765 

.921 

20 

.490 

.608 

.741 

.843 

.952 

.997 

1  .0 

50 

M.M.L. 

Method 

R  * 

85# 

Confidence 

Level 

Sample  Size 

.50 

.60 

.70 

.80 

.90 

•  95 

.99 

.300 

.300 

.300 

.300 

.354 

.448 

.635 

5 

.393 

.448 

.450 

.450 

.454 

.518 

.712 

10 

.287 

.386 

.533 

.647 

.647 

.647 

.739 

20 

.305 

.496 

.604 

.745 

.845 

.852 

.918 

50 
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V  Assessment  of  Method 


Results 

The  three- parameter  Double  Monte  Carlo  Method  was  run  for  1000 
repetitious  at  system  reliabilities  of  74%,  85%,  and  96%  with  component 
failure  sample  sizes  of  five,  ten,  twenty,  and  fifty.  Initially, 
seventy-five  points  were  used  for  the  component  sample  distributions  of 
reliability  estimates  to  be  compatible  with  previous  work  (Ref  15)»  The 
results  of  these  runs  are  listed  in  Table  X.  The  1000  repetitions  of 
each  combination  were  divided  into  six  runs  of  150  and  one  run  of  100  to 
keep  all  required  CPU  times  below  4000  seconds. 

In  order  to  check  the  sensitivity  of  the  Double  Monte  Carlo  method 
to  the  number  of  points  in  the  component  sample  distributions  of 
reliability  estimates,  the  number  of  points  was  increased  to  150  and  the 
method  used  for  a  system  reliability  of  85%  with  component  failure  sample 
sizes  of  five,  ten,  and  twenty.  The  results  of  these  runs  are  listed  in 
Table  XI. 

Since  increasing  the  number  of  points  in  the  component  sample 
distributions  of  reliability  estimates  did  not  increase  the  accuracy  of 
the  method,  the  number  of  points  was  decreased  to  fifty  and  the  Double 
Monte  Carlo  Method  was  checked  at  reliabilities  of  74%,  85%,  and  96%  with 
component  failure  sample  sizes  of  five,  ten,  and  twenty.  The  reduction 
in  the  number  of  points  from  seventy-five  to  fifty  resulted  in  a  1/3 
reduction  in  the  number  of  parameter  estimations  required.  Since 
parameter  estimation  takes  nearly  all  of  the  computer  CPU  time  for  the 
method,  this  reduction  in  points  also  resulted  in  about  a  1/3  reduction 
in  CPU  time.  This,  combined  with  limiting  the  failure  sample  sizes  to 


45 


TABLE  X 


Doable  Monte  Carlo  Results:  Distribution  Size  75 

R  =  74* 

Confidence  Level  Sample  Size 


R  =  85* 

Confidence  Level 


Sample  Size 


.50 

•  60 

.70 

.80 

.90 

•  95 

•  99 

.420 

.528 

.642 

.729 

.817 

.876 

.940 

5 

.470 

.578 

.675 

.757 

.848 

.910 

.976 

10 

.362 

.455 

.522 

.609 

.715 

.786 

.856 

20 

.522 

.628 

.730 

.842 

•  933 

.967 

•  989 

50 

R  = 

96* 

Confidence 

Level 

Sample  Size 

.50 

.60 

.70 

.80 

.90 

.95 

.99 

.371 

.501 

.607 

.714 

.807 

.872 

.934 

5 

.578 

.687 

.786 

.863 

.935 

.959 

.989 

10 

.111 

.150 

.201 

.290 

.445 

.570 

.763 

20 

•  566 

.  660 

.767 

.862 

.948 

.973 

.992 

50 

TABLE  XI 

Double  Monte  Carlo  Results:  Distribution  Size  150 


R  -  85* 

Confidence  Level 


Sample  Size 


.50 

.60 

.70 

.80 

.90 

.95 

.99 

.522 

.623 

.730 

.824 

.879 

.933 

5 

RV-'Ji 

.580 

.690 

.767 

.856 

.900 

.958 

10 

.405 

.486 

.570 

.663 

.743 

.802 

.880 

20 

20,  allowed  an  increase  to  500  repetitions  within  4000  seconds  of  CPU 
time.  The  results  of  these  runs  are  listed  in  Table  XII. 

The  results  indicate  that  a  sample  size  of  twenty  had  the  greatest 
bias  of  the  four  sample  sizes  tested.  This  was  particularly  noticeable 
at  a  system  reliability  of  96^.  To  obtain  a  better  assessment  of  the 
bias  trend,  the  Double  Monte  Carlo  method  was  run  for  a  system 
reliability  of  9b%  at  sample  sizes  of  fifteen,  twenty-five,  and  thirty. 
Fifty  points  were  used  for  the  sample  distributions  of  reliabilities. 

The  results  of  these  runs  are  combined  with  the  results  of  the  runs  under 
the  same  conditions  except  with  failure  sample  sizes  of  five,/  ten,  and 
twenty  in  Table  XIII  to  show  the  trend  in  the  system  reliability  bias. 

Discussion  of  Results 

Parameter  Estimation.  The  accuracy  of  the  Double  Monte  Carlo  method 
is  dependent  on  the  accuracy  of  parameter  estimation.  In  order  to  be 
practical  for  computerized  operation,  the  parameter  estimation  method 
used  must  also  be  fast.  If  the  bias  of  the  reliability  estimates  is 
known,  these  estimates  can  be  corrected  for  bias  in  the  program. 

Depending  on  the  parameter  estimation  technique  used,  the  bias  is  not 
necessarily  directly  related  to  sample  size.  In  the  case  of  the  L.L.S. 
method,  the  routine  is  fast  and  reasonably  accurate  but  has  no  particular 
optimum  properties;  therefore,  increasing  sample  size  does  not  imply  that 
the  bias  will  decrease.  The  results  indicate  that  the  system  bias 
becomes  increasingly  negative  as  the  sample  size  increases  from  five  to 
twenty,  then  increases  to  slightly  positive  for  a  sample  size  of  thirty 
(Table  XIII).  Generally,  the  bias  will  not  be  known  and  it  will  prove 
more  practical  to  either  establish  a  system  bias  or  use  the  confidence 
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TABLE  XII 


Double  Monte  Carlo  Results:  Distribution  Size  50 


R  »  85* 

Confidence  Level  Sample  Size 


.50 

.60 

.70 

.80 

.90 

.95 

.99 

5 

10 

20 

|___ _ — 

.452 

.464 

.596 

.352 

.561 

.481 

.627 

.655 

.556 

.721 

.758 

.659 

.817 

.855 

.757 

.877 

.915 

.816 

.955 
•  958 
.890 

R  = 

96$ 

Confidence 

Level 

Sample  Size 

.50 

.60 

.70 

.80 

.90 

.95 

.99 

mm 

n 

.717 

.806 

.874 

•  924 

5 

M.V-B2 

Efim 

.754 

.860 

.914 

.958 

10 

.138 

.186 

.241 

.555 

.482 

.622 

.785 

20 

R  = 

74$ 

Confidence 

Level 

Sample  Size 

.50 

.60 

.70 

.80 

.90 

.95 

.99 

.475 

.570 

.648 

.758 

.825 

.882 

.952 

5 

.472 

.565 

.664 

.755 

.856 

.912 

.957 

10 

.461 

.556 

.628 

.725 

.831 

.901 

•  942 

20 

TABLE 

XIII 

Extended 

Double 

Monte 

Carlo  Results: 

Distribution  Size  50 

R  = 

96$ 

Confidence 

Level 

Sample  Size 

.50 

.60 

.70 

.80 

.90 

.95 

.99 

.425 

.514 

.611 

.717 

.806 

.874 

.924 

5 

.455 

.571 

.654 

.754 

.860 

.914 

.958 

10 

.284 

.566 

.460 

.558 

.669 

.761 

.851 

15 

.150 

.186 

.241 

.555 

.482 

.622 

.785 

20 

.558 

.669 

.749 

.855 

.928 

.964 

.988 

25 

.554 

.655 

.755 

.822 

.925 

.962 

.994 

30 
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level  calculated  that  is  closest  to  the  desired  level.  Since  with  real 
data  the  exact  actual  parameters  will  be  unknown,  the  bias  or  levels  used 
must  be  for  parameters  close  to  the  true  parameters.  This  bias  can  be 
calculated  in  the  following  manner: 

1.  From  the  available  failure  data,  estimate  the  component 
parameters  as  accurately  as  possible.  Computer  CPU  time  is  not  a 
factor  for  this  estimation. 

2.  Use  these  accurate  estimates  as  true  parameters  and  run  a  Double 
Monte  Carlo  program  to  determine  accuracy  at  the  significance 
level(s)  desired.  Empirically  determine  system  bias  by  applying 
bias  to  the  system  reliability  estimates  until  the  results  are 
accurate.  This  bias  will  be  accurate  for  the  estimated  parameters 
(step  1 )  and  conditions  and  should  be  close  for  the  true  parameters. 
3-  Use  the  bias  (from  step  2)  to  calculate  the  appropriate 
confidence  limit(s). 

Sample  Distribution  Size.  Three  sizes  of  sample  distributions  of 
component  reliability  estimates  (50,  75,  150)  were  tested  with  no 
apparent  difference  in  accuracy  of  results.  The  actual  error  caused  by 
using  a  small  number  of  points  cannot  be  directly  calculated  since  the 
distribution  of  the  reliability  estimates  is  not  known  exactly,  but  since 
the  error  is  the  difference  between  the  point  on  the  true  curve  and  the 
point  obtained  by  interpolation  between  an  upper  and  lower  point  on  the 
true  curve,  the  error  will  be  much  less  than  the  spacing  of  the  points 
used  for  establishing  the  reliability  distribution.  The  use  of  50  points 
gave  good  results  in  these  trials  and  should  provide  the  required 
accuracy  for  any  application. 

General  Method.  The  Double  Monte  Carlo  method  of  confidence  limit 
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estimation  is  easy  to  implement  and  provides  good  consistent  results  for 
the  three- parameter  Weibull.  There  is  a  requirement  to  remove  bias  for 
maximum  accuracy,  but  an  estimate  of  bias  is  readily  obtainable  and  easy 
to  apply.  If  the  extra  accuracy  is  not  required,  then  considerable 
computer  time  can  be  saved  by  using  the  biased  estimators  for  approximate 
confidence  limits.  Four  parameter  estimation  methods  are  included  in 
Appendix  B,  but  the  method  can  be  readily  used  with  any  parameter 
estimation  technique.  In  the  results  of  this  development,  the  estimation 
of  the  location  parameter  appeared  to  be  the  most  critical;  therefore, 
the  data  should  be  examined  for  a  gross  estimate  of  this  parameter  before 
being  input  into  any  particular  routine.  A  good  approach  would  be  to 
ensure  the  data  is  well  spread  by  removing  some  portion  of  the  first 
order  statistic;  then  use  the  L.L.S.  method  of  parameter  estimation. 


i 

I 


50 


VI  Illustration  of  Method 


In  order  to  illustrate  the  method,  the  component  parameters  used  for 
the  system  reliability  of  0.96  were  used  to  generate  one  set  of  fifteen 
data  points  for  each  component  using  a  seed  value  of  17943*  These  data 
are  included  in  Appendix  D.  The  L.L.S.  method  was  used  to  estimate  the 
parameters  from  these  data  and  a  Double  Monte  Carlo  program  was  run, 
using  these  estimated  parameters  as  true  parameters,  to  calculate  system 
bias.  As  a  matter  of  interest,  the  Harter-Moore  three-parameter  maximum 
likelihood  routine  was  also  run  to  estimate  the  parameters  from  these 
data  points.  The  results  of  the  parameter  estimation  routines  are  listed 
in  Table  XIV. 

A  system  bias  of  -0.005  was  found  to  provide  good  results  at  the 
high  confidence  limits  of  interest  using  seed  values  of  7539  and  96  for 
the  random  number  generation  on  two  Double  Monte  Carlo  runs  of  500  each 
for  1000  total  runs.  As  a  test  of  the  accuract  of  this  bias,  it  was 
applied  to  the  system  reliability  of  1000  Double  Monte  Carlo  runs  using 
the  known  true  parameters  from  which  the  failure  data  were  generated  and 
seed  values  of  135  and  17*  Table  XV  lists  the  results  of  these  two  runs 
and  the  original  results  not  corrected  for  bias. 
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TABLE  XIV 


Parameters  for  Demonstration 


Component 

Method 

c 

k 

9 

1 

L  •  L  •  S  * 

69.5 

2.42 

226 

True 

10 

2.9 

270 

H-M  M.L. 

118 

1 .70 

173 

2  1 

L  •  L . S  • 

56.5 

1  .49 

458 

True 

0 

2.1 

470 

H-M  M.L. 

82.2 

1  .70 

400 

3 

L  •  L « S  * 

82.5 

2.36 

333 

True 

15 

2.3 

400 

H-M  M.L. 

.0000027 

3.63 

422 

4 

L  •  L .  S . 

46 

136 

True 

30 

1 

160 

H-M  M.L. 

.0000023 

183 

5 

L .  L .  S  . 

~  i 

80 

2.80 

206 

True 

25 

2.7 

250 

H-M  M.L. 

152 

1.59 

123 

6 

L .  L .  S . 

74.9 

1  .40 

229 

True 

50 

2.0 

280 

H-M  M.L. 

88 

1 .19 

213 

TABLE  XV 


Illustration  Results 


Confidence  Level 


.50 

.60 

.70 

.80 

.90 

.95 

•  99 

Bias 

Estimates 

.615 

.704 

.781 

.849 

.920 

•  949 

•  974 

True 

Parameters 

with  Bias  .335 

.413 

.507 

.604 

.711 

.795 

.866 

True 

Parameters 

without  Bias  .284 

.366 

.460 

.558 

.669 

.761 

.851 
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VII  Concluding  Material 

Summary  I 

A  Double  Monte  Carlo  method  of  estimating  complex  system 
reliabilities  at  any  confidence  level  was  developed  based  on 
three-parameter  Weibull  component  failure  distributions.  In  order  to 
implement  the  method,  three  fast,  accurate  parameter  estimation  routines 
were  developed  and  tested.  The  most  effective  routine,  the  L.L.S.,  was 
selected  and  used  in  the  Double  Monte  Carlo  program.  Good  results  were 
obtained  with  componet  failure  sample  sizes  as  low  as  five  for 
reliabilities  from  14%  to  96%.  A  step-by-step  procedure  with  an 
illustration  of  the  method  are  provided. 

Conclusions 

It  is  concluded  that: 

1.  the  Double  Monte  Carlo  method  can  be  readily  used  to  provide 
reliability  confidence  limits  for  complex  systems  with 
three-parameter  Weibull  failure  distributions; 

2.  sample  sizes  as  small  as  five  can  be  used  for  system  reliability 
confidence  limit  estimation; 

3-  the  Double  Monte  Carlo  method  is  not  adversely  affected  by 
reducing  the  number  of  points  in  the  sample  distributions  of 
reliabilities  to  fifty;  and 

4.  further  development,  full  testing,  and  establishing  biases  for 
three-parameter  Weibull  reliability  estimations  is  needed. 


* 
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Recommendations 


It  is  recommended  that; 

1 .  further  development  of  computerized  reliability  estimation 
routines  be  done;  and 

2.  bias  tables  be  developed  for  a  fast,  accurate  three-parameter 
Weibull  reliability  estimation  routine. 
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APPENDIX  A.  Double  Monte  Carlo  Program 


PROGRAM  SYSREL 

C  THIS  PROGRAM  PROVIDES  CONFIDENCE  LIMITS  ON  THE 
C  RELIABILITY  OF  A  COMPLEX  SYSTEM  WHOSE  COMPONENT  RELIABILITIES 
C  CAN  BE  ESTIMATED.  IT  ASSUMES  THAT  THE  UNDERLYING 
C  COMPONENT  DISTRIBUTIONS  ARE  3-PARAMETER  UEI BULL 
C  AND  USES  THE  DOUBLE  MONTE-CARLO  TECHNIQUE  TO 
C  GENERATE  RELIABILITY  ESTIMATES.  INTERNATIONAL 
C  MATHEMATICS  AND  STATISTICS  LIBRARY  CIMSLI  ROUTINES 
C  GGUBSt FOR  GENERATING  A  UNIFORM  (0*11  RANOOM  VARIABLE* 

C  AND  GSUIB*  FOR  GENERATING  A  ONE-PARAMETER  UEIBULL 
C  RANDOM  VARIABLE*  ARE  USED.  A  LEAST  SQUARES  METHOD 
C  OF  ESTIMATING  THE  THREE 

C  UNKNOWN  PARAMETERS  IS  USED  IN  ROUTINE  PARAM. 

C  THE  PROGRAM  GENERATES  ITS  CUN  REAL  FAILURE  DATA 
C  BUT  CAN  READILY  BE  MODIFIED  TO  USE  ACTUAL  FAILURE 
C  DATA  BY  DELETING  THE  APPROPRIATE  SECTION. 

C  inputs: 

C  REPS  -  THE  NUMBER  OF  TIMES  SIMULATED  RELIABILITIES  ARE  GENERATED 

C  TO  ESTABLISH  A  RELIABILITY  DISTRIBUTION 

C  FUNCTION. 

C  CCI> *KCI>,THETA(I>=  THE  PARAMETERS  OF  THE  COMPONENT 
C  FAILURE  DISTRIBUTIONS. 

C  TNFA I L=  THE  NUMBER  OF  FAILURES  OF  EACH  COMPONENT 

C  ........... 

c  outputs: 

C  FOR  50*60*70,30.50*95.99*  THE  NUMBER 

C  OF  TIMES  THE  LOUER  LIMIT  WAS  LESS  THAN 

C  THE  LOUER  CONFIDENCE  LIMIT  CALCULATED. 

C  *********** 

c  variables: 
c  integer: 

C  REPS=  AS  ABOVE 

C  TNFA  IL  =  AS  ABOVE 

C  NLL50=NUMBER  OF  TIMES  THE  50  PERCENT  LOWER  LIMIT 

C  CAPTURES  THE  RELIABILITY. 

C  NLL60  »  NLL7  0.NLL80* NLL9U «NLL95«NLL99-  AS  NLL50 

C  FOR  THE  HIGHER  LIMITS 

c  real: 

C  C=  DISPLACEMENT  PARAMETER 

C  K=  SHAPE  PARAMETER 

C  TiET A=  SCALE  PARAMETER 

C  TREL=  TRUE  RELIABILITY  OF  THE  SYSTEM 

C  TFAILI  TO  TFAIL6=  TRUE  FAILURE  VECTORS  FOR  THE 

C  COMPONENTS 

C  FCCn.FK(I).FTHETACI)=  FIRST  ESTIMATED  PARAMETERS 

C  OF  COMPONENT  I 

C  SC<I  >,SKCI>,STHETA(I>  =  AS  ABOVE  BUT  S ECOND  ESTIMATE 

C  SFAIL1  TO  3FAIL6=  SIMULATED  COMPONENT  FAILURES 

C  R  A  NX  a  >  =MEDI  AN  RANK 

C  RU=  VECTOR  OF  UNIFORM  RANOOM  NUMBERS 

C  SRI  TO  SR6=  SIMULATED  RELIABILITY  OF  EACH  COMPONENT 

C  SR  STIMULATED  RELIABILITY  OF  SYSTEM 
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FIRSTS  EXTRAPOLATED  VALUE  CORRESPONDING  TO  THE 
FIRST  MEDIAN  RANK  OF  0 
LAST=  AS  ABOVE  BUT  LAST  MEDIAN  RANK  OF  1 
SR C=  ORDERED  RELIABILITY  OF  EACH  COMPONENT  WITH 
EXTRAPOLATED  END  VALUES 
SASORD=AS  SRC  BUT  FOR  SYSTEM 
CL 50 . CL99s  LOWER  CONFIDENCE  LIMITS 


INTEGER  REPS*TNFAILfNLL50,NLLSO*NLL70.NLL90 
INTEGER  NLL90  * NLL95  * NLL99 
INTEGER  RGEUNRUNS 

REAL  R(6)*C(6)*K(&  )*THETA(G)*TREL*  FC(6)»FK(5) 

REAL  FTHEIAC6) 

REAL  TFAIL1 C 1 00 ) » TF AIL2 ( 100 ) *  TF  A IL3 < 10 0 ) *TFAIL A C 1 0 0 ) * TF AI L5 ( 100  ) 
REAL  TFAIL6C100)*SFAIL1(100 ),SFAIL2C100)*SFAIL3C100)*$FAILA(100) 
REAL  SFAIL5C100 )*SFAIL6(100 ) *SC(6)*SK(6) *STHETACS)»RANK(602) 

REAL  SRS0RDC602 ) *  P  U  (  6  ) 

REAL  FIRST.LAST.SRSC600 ), SRI (300) *SR2C300) 

REAL  SR 3(300) *SR4 C3 00) *SR5C 300) *SR6 C300 ) *SRC(£*332) 

REAL  CL50*CL60,CL70 *CL80,CL90 
REAL  CL95  *CL99 
DOUBLE  PRECISION  DSEED 
DECLARATIONS  COMPLETE 


INITIALIZE  INPUTS  AND  COUNTERS  FOR  NUMBER  OF  TIMES 
LOWER  LIMIT  EXCEEDED 
DSEED=179A3.0D0 

PRI NT***L.L.S.  DSEED  =  **OSEEO 
RGElsO 
REPS= 150 
N  RUNS =150 

PRINT***  NRUNS  =  »»NRUNS 

TNFAIL=10 

CCI  )  =  10. 

CC2  )s0. 

C (3 )= 15* 

CCA )=30. 

C (5 ) =25» 

CCS )s50. 

K (1 )=2.3 
K  (2  )  =  •  8 
K (3 )=2.0 
K(A)s3.5 
K  (5  )  =  1  ■  2 
KCS)=2. 

THETAC1 )=200. 

THET AC  2 ) =A70« 

THETAC3)sl80. 

THETAC A)s120. 

T  HET  A (5 ) =250 • 

THCTAC6)=150. 
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N  LL50  =0 
NLL60-0 
NLL70=0 
NLL80=0 
NLL*0  =  0 
M  LL.  85=0 
NLL9*=0 

CALCULATE  THE  POINT  ESTIMATE  OF  TRUE  SYSTEM  RELIABILITY 

CALL  RLBTYCC*K»THETA*TREL*RGE1I 
C  GENERATE  THE  REQUIRED  NUMBER  OF  MONTE  CARLO 
C  REPETITIONS.  FOR  EACH  RUN  DETERMINE  IF  THE  TRUE 
C  RELIABILITY  IS  ABOVE  THE  LOWER  LIMIT  AT  EACH 
'  C  CONFIDENCE  LEVEL  AND  INCREMENT  THE  APPROPRIATE 
C  COUNTER  IF  SO. 

DO  5  H=1 » NRUNS 
C 

C  GENERATE  THE  TRUE  COMPONENT  FAILURE  DATA  FOR  EACH 
C  COMPONENT 

CALL  UEIBLCDSEED*C< 1 ) *K (1 >* T HET AC  1 >  * TNFA IL* TFA I  LI  I 
CALL  WEIBL(0SEED»C(2)  »K (2  >• THE T AC  2) • TNFA IL»T FAI L2> 
CALL  UEIBL(DSEED*C(3)*K(3)*THETA(3) *  TNFAIL*  TFAI L3 1 
CALL  UEIBL<DSEED*C< ♦  > »K (♦ >* THETA<4>  * TNFA IL* TFAI L4> 
CALL  VEIBLCDSEED.CC  5)  »K <5 )* THET A< 5) •TNFAIL* TFA I LSI 
CALL  UEIBL(DSEED*CC6)*KC6)*THETA<6>  *TNFAIL*TFAILS> 

C  SORT  THE  TRUE  FAILURE  DATA  AND  CALCULATE  THE 
C  ESTIMATORS  OF  THE  PARAMETERS 
C 

CALL  VS»TACTFAIL1.TNFAIL> 

CALL  VSRTACTFAIL2.TNFAIL) 

CALL  VSRTA(TFAIL3«TNFAIL) 

CALL  VSRTACTFAILA *TNFAIL> 

CALL  VSRTAC TFAIL5fTNFAIL> 

CALL  VSRTAtTFAIL6*TNFAIL) 

CALL  PARA  M(TNFAIL*TFAIL1«FC(1)*FK(1)*FTHETA<1>) 

CALL  PA RAHC TNFAIL  *TFAIL2*FCC2>*FK<2)» FT HET A<2>) 

CALL  PARA MC TNFA IL»T FA IL3*FC < 3) »FK C3 > *FTHETA< 3 >» 

CALL  PARA M( TNFA I L*T FA  IL4  *FC (4)»FK(4)*FTHETA(4)I 
CALL  PARAM«TNFAIL»TFAIL5»FCC5>»FKC5)*FTHETAf5)> 

CALL  PARA  M(TNFAIL*TFAIL6*FC(6)*FK(6) * FT HET  AC  S I  > 

C  FOR  THE  EMPIRICAL  DISTRIBUTION*  GENERATE  REPS 
C  RELIABILITY  ESTIMATES. 

C  FIRST  GENERATE  TNFAIL  SIMULATED  FAILURES  FOR  EACH 
C  OF  THE  COMPONENTS  USING  THE 
C  ESTIMATES  OF  THE  PARAMETERS 
DO  10  L  =  l *REP S 

CALL  UEIBL( DSEED*FC (1 ) *FK (1 > »FTHETA C 1 1  * TNFAIL * SFAI LI  ) 
CALL  UEIBL(DSEED*FC(2)«FK(2)*FTHETA(2)*TNFAIL«SFAIL2) 
CALL  UEIBL(DSEED*FC(3)*FK(3)*FTHETA(3) «TNFAIL*SFAIL3> 
CALL  UEIBL(DSEED*FC(A)*FK(4)*FTHETA(4)*TNFAIL*SFAILA) 
CALL  UEIBL(DSEED*FC<5)*FK<5>*FTHETAC51*TNFAIL*SFAIL5» 
CALL  WE IBL( DSEED* FC  (6 ) *FK (6 ) *FTHETA (6 ) • TNFAIL *SFA IL6> 
C  ORDER  THE  SIMULATED  FAILURES  AND  USE  THEM  TO  OBTAIN 
C  THE  SECONO  ESTIMATE  OF  THE  COMPONENT  PARAMETERS 
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CALL  VSRTACSFAIL1 *TNFAIL) 

CALL  VSRTACSFAIL2»TNFAIL) 

CALL  VSRTACSFAIL3»TNFAIL> 

CALL  VSRTAC$FAIL4,TNFAIL> 

CALL  VSRTACSFAIL5,TNFAIL1 
CALL  VSRTACSFAIL6*TNFAIL> 

CALL  PARANCTNFAIL*SFAIL1»SCC1)*SKCI),STHETACI>> 

CALL  PARA MC  TNFA IL*SFA IL2»SC C2)*SKC2)»STHETAC2>> 

CALL  PARAMCTNFAIL*SFAIL3*SCC3>* SKC3 > * STHETAC 3 > > 

CALL  PARAMCTNFAIL*SFAIL4*SCC4)*$KC4)*STHETAC4)) 

CALL  PARAHC TNFA IL*S FAILS. SC C 5) *SKC5) • STHETAC5) > 

CALL  PARAMCTNFAIL»SFAIL5»SCC6)*SKC6I*STHETAC6)> 

C  CALCULATE  THE  RELIABILITY  OF  EACH  COMPONENT 
C  USING  THE  SECOND  ESTIMATE  OF  COMPONENT  PARAMETERS 
C  TO  BUILD  THE  COMPONENT  RELIABILITY  VECTORS  SRI  TO  SRS 
C 

SRI  (L)=RELYCSCC1>*SK(1)« STHETAC 1) *RGE1) 

SR2 (L  )=RELY CSCC2)tSKC2)*STHETAC2)»RGEl> 
SR3(L)=RELY(SC(3) »SKC3)»STHETAC3) «RGE1 ) 
SR4(L)=RELYCSCC4)*SKC4)*STHETAC4) *R GE1 > 
SR5CL)=RELYCSCC5>  «SKC5) , STHETAC 5) »RGEl ) 
SRGCL)=RELYCSCC6>*SKC6)*STHETAC6>«RGE1 > 

10  CONTINUE 
C 

C  ESTABLISH  RELIABILITY  DISTRIBUTION  FUNCTIONS  FOR 
C  EACH  COMPONENT  TYPE  USING  ORDERED  MEDIAN  RANKS  I 
C  ONE  RANK  VECTOR  AND  THE  SIX  RELIABILITY  DISTRIBUTION 
C  VECTORS  MILL  MAKE  UP  THE  SIX  DISTRIBUTION  FUNCTIONS 
C  EACH  VI TH  REPS>2  VALUES 
C 

RANK!  1)=0  . 

R  ANM  REPS *2  )  =  1« 

DO  15  I=i tREPS 

RANKCI*l)=CREALCI)-.3)/CREALCREPS)*.4) 

CONTINUE 

ORDER  THE  RELIABILITY  VECTORS  AND  ESTABLISH  THE  VALJiS 
CORRESPONDING  TO  MEDIAN  RANKS  0  AND  1 

CALL  VSRTACSR1.REPS  > 

CALL  VSRTACSR2.REPS > 

CALL  V3RTACSR3.REPS > 

CALL  VSRTACSR4.REPS I 
CALL  VSRTA<SR5*REPS  I 
CALL  VSRTACSR6*REPS> 

CALL  EX  TRAC  SRI* RANK »FIRST*LAST*REPS> 

SRC (1 *1)=FIRST 
SRC  Cl  *REPS«-2)=LAST 

CALL  EX  TR  AC  SR 2*  RANK  *FIRST*L  AST*REPS ) 

SRCC2*1>=FIRST 
SRCC2*REPS*2>=LAST 

CALL  EXTRACSR3*RANK  *FIRST«L AST* REPS) 

SRC C3*l>= FIRST 
SRC  C3»REP3»2)=LAST 

CALL  EXTR ACSR4»RANK *F1RST«LAST*REPS) 

SRC (4*1 )=FIRST 
SRC  C4*REPS*2)=LAST 

CALL  EXTRACSR5*RANK *FIRST*LAST*REPS> 

SRC  C5  *1 )  =  FIPST 
SRCC5*REPS*2)=LAST 
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CALL  EXTRACSR6 »R  ANK  ,FIRST,L AST, REPS > 

SRC  C6,l>=FIRST 
SRC  C6,REPS*2>=LAST 
00  20  1=1 f REPS 
SRC  Cl ,1*1  )  =  SR1 ( I ) 

SRC  C2,I»1  )=SR  2< I > 

SRCC3*I*1)=SR3CI) 

SRC  C  4,IM  >  =  SR4C  I I 
SAC  C5,I>1 >  =  SR5CI I 
SRC (6«I*1  >=SR6CI> 

20  CONTINUE 

C  RANDOMLY  SELECT  A  RELIABILITY  FROM  EACH  RELIABILITY 
C  DISTRIBUTION  FUNCTION  USING  A  UNIFORMCO,l>  GENERATOR 
C  TO  OBTAIN  A  MEDIAN  RANK*  AND  LINEAR  INTERPOLATION. 

C  CALCULATE  SIMULATED  SYSTEM  RELIABILITY  SRS.  REPEAT 
C  600  TIMES. 

C 

00  25  1  =  1 ,600 

CALL  GGUBSCDSEED*6,RU> 

00  30  J=1 ,6 
11  =  2 

35  IFC II.LE.REPS*2»THEN 

IFCRUC JJ.LE.RANKCII) >THEN 

CALL  INTERPCRANKCII>,RANKCII-1>,SRCCJ,II>*SRCC J.II-1 ), 
IRUC  J)  *R  (J  )) 

1 1 =REPS*3 
ELSE 
II  =  II«-1 
END  IF 
GO  TO  35 
END  IF 

30  CONTINUE 

SRS  Cl  >=RC1>*CCCRC2I«-RC3>-RC2>*RC3>>*AC4>>«-CRC5>«-|*C6> 

1-AC  5)*RC6>)-C CRC2>fRC3>-RC2>*RC3>)*RC4>>* 

1CRC  5I^RC6»-R<5)  *R  C6II1 
25  CONTINUE 

C 

C  ORDER  THE  SYSTEM  C600>  RELIABILITY  ESTIMATES  USING 
C  MEDIAN  RANKS  AND  DETERMINE  THE  99,95*90*80,70,60,50 
C  PERCENT  LOWER  LIMITS.  NOTE  IF  EACH  CONTAINS  THE  TRUE 
C  SYSTEM  RELIABILITY  AND  IF  SO*  INCREMENT  THE  APPROPRIATE 
C  COUNTER  NLL50*...NLL99. 

CALL  VSRTACSRS»600) 

R  AMKC  1)  =  0  • 

A  Alt  K(  6021=1 . 

00  40  1=1*600 

RAMK(I+1)=(REAL<I >-.31/600.4 
40  CONTINUE 

CALL  EX TRAC SRS, RANK ,FIRST,LAST ,600 > 

SRSORDC 1>=FIRST 
SRSORDC602) =LAST 
00  45  1=1,600 
SRS0R0CI«-1>=SRSCI  > 

45  CONTINUE 


63 


CALL  INTERP  <R ANKC  8)  «R  ANKC  7)  *  SRSOR  0(8  ) *SRSORD<  7) « 
1.01.CL99) 

CALL  INTERP(RANK(32)*RANK(31)*SRS0RD(32)*SRS0RD(31)* 
1.05.CL95) 

CALL  INTERP (RAN K( 62  )  .RANK  <61 ) .SRSORD (62) *SR30RD(& 1) • 
1.1* CL90) 

CALL  INTERP(RANK(122)*RANK< 121) *  SRS0RDU22)  «SR  SORO<  121  > 

1.2, CL30> 

CALL  INTERP(RANK(182>*RANK< 181) tSRSORO< 182)* SRSOR 0<181 > 

1.3. CL70) 

CALL  INTERP <RAMK< 242 >  *RANK< 241) *SRS0RD(242>* SRSOR 0(241) 
1.4*  CL60) 

CALL  INTERP(RANK(302)*R  ANKC  301 ) • SRSORD (302). SR SORO( 301 ) 
1.5.CL50) 

I F< CL99.LT.TREL)NLL99=NLL99*1 
IF(CL95.LT.TREL)NLL95=N LL95+1 
IF(CL90.LT.TREL)NLL90=NLL90*i 
I F< CL80.LT.TREL)NLL30=NLL80*1 
IF(CL70.LT.TREL)NLL70  =  NLL70*1 
IF(CL60.LT.TREL)NLL60=NLL60*l 
IF< CL50.LT.TREL)NLL50=NLL50*1 
5  CONTINUE 

C  PRINT  THE  NUMBER  OF  TIMES  THE  TRUE  RELIABILITY 
C  UAS  ONER  EACH  CONFIDENCE  LIMIT 

PRINT***  TRUE  RELIABILITY  IS  •  • TREL 
PRI  NT** 'NUMBER  OF  FAILURES  *  * TNFAI L 
PRINT***  NUMBER  OF  REPETITIONS  *»REPS 
PRINT** 'NUMBER  OF  TIMES  RELIABILITY  GE  t  =  **R6E1 

PRINT***  NUMBER  ABOVE  50  PERCENT  LOWER  LIMIT  **NLL50 

PRINT*. •  NUMBER  ABOVE  60  PERCENT  LOWER  LIMIT  *,NLL60 

PRINT***  NUMBER  ABOVE  70  PERCENT  LOWER  LIMIT  **NLL70 

PRINT***  NUMBER  ABOVE  80  PERCENT  LOWER  LIMIT  *«NLL80 

PRINT***  NUMBER  A80VE  90  PERCENT  LOWER  LIMIT  **NLL90 

PRINT***  NUMBER  ABOVE  55  PERCENT  LOWER  LIMIT  **NLL95 

PRINT***  NUMBER  ABOVE  99  PERCENT  LOWER  LIMIT  **NLL99 

STOP 
END 

. . **** 

SUBROUTINE  EXTR A< X* Y .FI RST* LAST *REPS ) 

C  USES  L I NEAR  INTERPOLATION  OFF  THE  TWO  ENO 
C  VALUES  TO  OBTAIN  THE  RELIABILITIES  CORRESPONDING 
C  TO  THU  MEDIAN  RANKS  0  AND  1*  FIRST  AND  LAST. 

C  SLOPE  IS  THE  SLOPE  OF  THE  LINEAR  LINE  BETWEEN  THE 
C  TWO  END  VALUES  AT  EACH  END. 

INTEGER  REPS 

REAL  X(300)*Y(302)*FIRST»LAST 
Z  =X  (  2  )—X(l  ) 

ZZ=  X  ( REPS )-X< REPS-1  ) 

V=Y ( 3  )-Y<  2 ) 

VV=Y(REPS*1 )-V<  REPS ) 

I F( Z.GT • 0  .) THEN 
SLOPE= V/Z 

IF(SLOPE.EQ.O.)THEN 
FIR  ST  =  X ( 1 ) 


ELSE 

FIRST=XC1)-Y(2) /SLOPE 
END  IF 
ELSE 

FIRST=X<1I 
END  IF 

IFf  Z2. GT.O.JTHEN 
SLOPE= VV/ZZ 
IFCSLOPE.EQ.O. >T HEN 
LAST=XCREPS> 

ELSE 

LAST=X<REPS)*C<l-Y<REPS*l> > /SLOPE) 

END  IF 
ELSE 

LAST  =X (REPS) 

END  IF 

I FC  FIRST .LT.0.)FIRST=0. 

IFCLAST.GT.l. )LAST=1. 

RETURN 

END 

(;•*•*  ********************  *********  ****** 

SUBROUTINE  I NTERPC UPPER Y* LOWERY vUPPERR (LOWER R fHEDtR > 
C  USES  LINEAR  I NTERPOLAT I  ON  TO  OBTAIN  A  RELIABILITY 
C  Rt  CORRESPONDING  TO  THE  MEDIAN  RANK t  RED.  MEO  IS 
C  BRACKETED  BY  MEDIAN  RANKS  UPPERY  AND  LOWERY  WHICH 
C  CORRESPOND  TO  RELIABILITIES  UPPERR  AND  LOUERR. 

REAL  UP  PE  RY  (LOWER  Y»  UPPERR (LOUERR  t  R  (SLOPE  t  MED 
X  =UPPERR-LOUERR 
Y=UPPERY-LOUERY 
I F( X.GT  *0 •) THEN 
SLOPE= Y/X 

R=LOUERR*((MED-LOWERY>/SLOPEI 

ELSE 

R= UPPERR 
END  IF 

IF( R.GT.l.)RsUPpERR 

RETURN 

END 

C ******  ************** **•••**•*** 

SUBROUTINE  UEIBLCDSEED*C(K(THETA(NU(RW> 

C  CALCU.ATES  RANDOM  3-PARAMETER  UEIBULL  VARIATES 
C  RU.  A  TOTAL  OF  NW  ARE  FOUND  USING  IMSL  ROUTINE 
C  GGWIB 

C  C  IS  THE  POSITION  PARAMETER(  K  IS  THE  SHAPE (  AND 
C  THETA  IS  THE  SCALE. 

INTEGER  NW 

REAL  RUf300)(C(K( THETA 
DOUBLE  PRECISION  DSEED 
CALL  GGWIB<OSEED(K(NW(RW> 

00  3  1=1 (NW 
RUC  I >=THETA*RW< I )*C 
3  CONTINUE 

RETURN 
END 
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C ******  ************ 

SUBROUTINE  RLBTYIC*  KtTHET A*RS *RGEl ) 

C  CALCULATES  SYSTEM  RELIABILITY  FROM  COMPONENT 
C  PARAMETER  DATA. 

C  RC6)  IS  THE  COMPONENT  RELIABILITIES?  RS  IS  THE  SYSTEM. 
REAL  CC6) *K <6 ) *THET A(6)«R(6)«RS 
INTEGER  RGE1 

C  CALCULATE  COMPONENT  RELIABILITIES 
DO  10  1=1*6 

R (I  >  =RELY (C(I)fKCI)* THETAC I ) *RGE1 ) 

10  CONTINUE 

C  CALCULATE  SYSTEM  RELIABILITY 

RS=R<i)*(((R(2>*RI3)-R(2)*R<3>l*RCAI)*CRC5)+R<6l 
l-R< 5 >*R <*>)-( (R(2)+R(3>-R(2)*R(3))*R(4>>* 

1(R(  S)*R(6>-ft(5)*R(6)>) 

RETURN 

END 

C****** ************* 

FUNCTION  RELY(C»K*T HET A *RGE 1  ) 

C  CALCULATES  COMPONENT  RELIABILITY  FOR  3-PARAMETER  UEI  BULL 
RE4L  C*K.THETA»T 
INTEGER  RGE1 
T=100. 

IF( X*LE. 0 .) THEN 
RELY  =  t • 

RGE1  =  R  GE 1*1 
ELSE 

X=(CT-C)/THETA)**K 

IFCX.LT.20.)THEN 

RELY=EXP<-X> 

ELSE 

REL Y=  0 • 

END  IF 
E»'D  IF 
RETURN 
END 
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APPENDIX  B.  Parameter  Estimation  Routines 


C************************** 

SUBROUTINE  PARAMCNUM*FAILtCU.KU* THETAU) 

C  MAXIMUM  SAMPLE  SIZE  3IMENSIONED=100 
C  INPUT 

C  NUM=SAMPLE  SIZE 
C  FAIL  =V ECTOR  OF  SAMPLE  DATA 
C  OUTPUTS 

C  ESTIMATES  OF  LOCATION  PARAMETER  CU*  SHAPE  PARAMETER 
C  KU,  AN D  SCALE  PARAMETER  THETAU 

c  variables: 

C  INTESER 
C  NUM=SAMPLE  SIZE 
C  MIN=MI  NIMUM  SQUARE  ERROR 

C  FLAG  IS  A  FLAG  TO  MARK  THE  MINIMUM  ERROR 
C  REAL 

C  FAIL=VECTOR  OF  SAMPLE  DATA 

C  RANK=VECTOR  OF  MEDIAN  RANKS 

C  PLOTRK  =PLOTT ING  RANK  =  LN( LNC l -1 /F< T ) » > 

C  EY  =  VICTOR  OF  EXPECTED  Y  VALUES  IF  THE  SAMPLES  PLOT 
C  LINEARLY  ON  A  UEIBULL  PLOT  WITH  THE  Y  AXIS  BEING 

C  THE  LN  OF  THE  FAILURE  TIMES  AND  THE  X  AXIS  THE 

C  LNLN  OF  1/C1-F(T)» 

C  SQRERR  =VECTOR  OF  SQUARE  ERROR  SUMS 
C  LNFA IL=VECT0R  OF  LN  OF  SAMPLE  DATA 
C  CU*KU* THETA W=  AS  ABOVE 

C  EC*EK»ETHETA=  ESTIMATES  OF  PARAMETERS 
C  T  0  T  A  LX  =  THE  TOTAL  OF  THE  ABSCISSA  VALJES 
C  TOTALY=  THE  TOTAL  OF  THE  ORDINATE  VALUES 

C** *••*.** •«•**•* **•**•*••**•«••* ******* 

C 

INTEGER  MUM,  FLAG 

REAL  FAIL  Cl  00 > .RANK <100 > *PLOTRK <100 )*EY<0:10> 

REAL  SQRERR(0:i0>«LNFAIL<100  >*CU*KU* THETAU* MIN 
REAL  EC(0:i0>  *EK(0: 10 ) tETHETACO 1 10) *TOTALX* TOTAL f 
C  DECLARATIONS  COMPLETE 
C********  ******************** 

C  CREATE  A  MEOIAN  RANK  VECTOR  AND  A  PLOTTING  LNLN 
C  VECTOR  WITH  A  TOTAL  FOR  THE  X  AXIS 
TOTALX=0. 

DO  10  1=1 »NUM 
R ANKCI)=<I-.3>/CNUM*.A) 

PLOTRKC I)=ALOGCALOG (1 •/ <1 .-RANK C I ) ) >> 
TOTALX=TOTALX*PLOTRKC I > 

10  CONTINUE 
C 

C  SET  VALUES  OF  C  FROM  0  TO  l.*FAIL(l)  AND 
C  DETERMINE  THE  SQUARE  ERROR  FROM  A  LINEAR 
C  PLOT  3 N  UEIBULL  PAPER  FO»  EACH 

C  USIN6  THE  RELATIONSHIP  LNC T) =C l/K) *LNLNC1/1-FCT >> 

C  *LN<  THETA) 


67 


-  Z  =PLOTRKC  NUMJ-PLOTR  KC 1 > 

90  15  J=0.10 
N=l 

IFC  J.EQ.10JTHEN 
N= 2 

2=PL0TRKCNUMJ-PL0TRKC2> 

ENO  IF 

ECCJ)=FAILC1)*.1*R£ALCJ) 

00  20  I =N» NUM 

LNFAILC IJ=AL0GCFAIL  Cl »-ECC J > ) 

20  CONTINUE 

C  FIND  THE  RELATIONSHIP  THAT  IS  CLOSEST  TO  LINEAR 
C  AND  USE  THE  PARAMETERS  FOR  THAT  MATCH 
C  FIRST  ESTIMATE  K*  THEN  THETA  ASSUMING  THAT 
C  1/K  IS  THE  SLOPE  OF  THE  LINE 

EKCJ)=Z/C LWFA IL( WUM  >-LNFAIL C  N ) ) 

TOT  ALT=0. 

DO  25  I=N*NUP 
T OT  AL Y= TOTAL Y*LNF AIL  C 1 1 
25  CONTINUE 

ETH  ET  AC  J) =CTOTALY-TOTALX/EK  C  J) > /PEAL  <NUM-N*1 ) 

C  CALC U. ATE  THE  ESTIMATES  OF  LNCT-CI  UHICH  IS  THE  ORDINATE. 

C  ESTIMATE  THE  ERROR  SQUARED 
SQRERRC  J)=0. 

DO  50  I'MvNUn 

ETC  I )=PLQTRKC IJ/EKC  JJ+ETRETAC J> 

SQRERRC JI=SQRERR(J) ♦(EYCI)-LNFAIL(I ) >  MEY Cl  J-LNF AILC I  > » 
30  CONTINUE 

15  CONTINUE 

C  FIND  THE  CLOSEST  LINEAR  RELATIONSHIP  ANO  USE  THE 
C  PARAMETERS  AS  THE  ESTIMATES 
FLAC  =  0 

M IN = SQRERRC 0) 

00  35  1=1*10 

IFC  SQRERRCIJ.LT. MIN) THEN 
HI  N=SQRERP  C I ) 

FLAG=I 
END  IF 

35  CONTINUE 

CW=ECCFLAG) 

KU=EKCFLAG) 

THETA«=EXPCETHETACFLAGI j 

RETURN 

ENO 
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c**«* •»•»**•»*****•*****•*• 

SUBROUTINE  PARAMCNUM*FAILtCU*KU» THETAU) 

C  MAXIMUM  SAMPLE  SIZE  DIMENSIGNED=100 
C  INPUT 

C  NUM=SA  MPLE  SIZE 
C  FA IL  =V ECTOR  OF  SAMPLE  DATA 
C  OUTPUTS 

C  ESTIMATES  OF  LOCATION  PARAMETER  CU»  SHAPE  PARAMETER 
C  KU»  AN  0  SCALE  PARAMETER  THETAU 
C  VARIA3LESI 
C  INTESER 
C  NUM=  SAMPLE  SIZE 
C  3  EAL 

C  FAIL=VECTOR  OF  SAMPLE  DATA 
C  RANK=VECTOR  OF  MEDIAN  RANKS 
C  PLOTR<=PLOTTING  RANK  =  LN< LN <1 -1/F <T ) ) ) 

C  LNFA IL  =  VECTOR  OF  LN  OF  SAMPLE  OATA 
C  CU»KU*  THETAW=  AS  ABOVE 

C  TOTALX  =  THE  TOTAL  OF  THE  A8SCISSA  VALUES 
C  TOTA  LY  =  THE  TOTAL  OF  THE  ORDINATE  VALUES 


C 

INTEGER  NUM 

REAL  FAIL  Cl  00) »RANK CIOQ > *PLOTRK< 1 03 ) 

REAL  SLOPE tLNFAILCl 01) tCU*KUt THETAUfTOTALXt TOTAL Y 
C  DECLARATIONS  COMPLETE 
C***»***» • •••*»•*•»••••»***•* 

C  CREATE  a  MEDIAN  RANK  VECTOR  AND  A  PLOTTING  LNLN 
C  VECTOR  WITH  A  TOTAL  FOR  THE  X  AXIS 
TOT ALX=0. 

DO  10  I=1*NUM 

RANK (I  )  =  CI-.3)/CNUM*.A) 

PLOTRICC  I)=AL0GCALCGC1./C1.-RANKCI))  )) 

TOT ALX=TOTALX»PLOTRKCI) 

10  CONTINUE 

C  ESTIMATE  C  BY  LINEAR  EXTRAPOLATION  FROM  THE  FIRST 
C  TWO  ORDER  STATISTICS 

CU=  F  A IL C2)-CFAILC2)-FAIL<1) )*RANK < 2  ) /< RANK <2 >-R  ANiCf  1 )  > 
C  CALCU. ATE  THE  LOG  OF  THE  FAILURE  TIMES  MINUS  C 
C  AND  THE  TOTAL  FOR  THE  Y  AXIS 
TOT ALY=0* 

DO  25  1=1 »NUM 
LNFAILCI)=ALOG(FAIL  <I)-CW) 

TOTALY=TOTALY«-LNFAILCI) 

CONTINUE 


25 


C  CALCULATE  K  BY  THE  AVERAGE  SLOPE  INVERTED 
C  DETERMINE  THE  SUM  CF  ALL  THE  SLOPES  AND  THE  AVERAGE 
C  INVERTED 

SLOPE=0. 

DO  20  N=2 »NUH 

S  LO  PE  =  SLOPE ♦( LNFAIL (Hl-LNFA IL  (M-l »l /< PLOTRK(H>-PLOTRK (M-l > > 
20  CONTINUE 

KU=  <R£AL(NUM>-1 •) /SLOPE 

C  LNCTHETA)  IS  A  CONSTANT  WHICH  IS  INCLUDED  IN  THE 
C  THE  RE_  ATIONSHI P  LMT  J=1/K*LNLN(1/<I-F< T)l > 

C  ♦  LN(THETA).  IN  THE  CALCULATIONS  OF  THE  TOTALS 
C  IT  IS  INCLUDED  NUM  TIMES*  SO  THE  AVERAGE  VALUE 
C  WILL  E3UAL  LN(THETA) 

THE  TAW- (TOTAL Y- TO TALX/KW) /REAL (NUR) 

THETA W~EXP<THETAU> 

RETURN 

END 
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C****** ************** ****** 

SUBROUTINE  PARAM<NUM*FAIL*CW*KW,THETAy) 

C  MAXIMUM  SAMPLE  SIZE  DIMENSIONED=100 
C  INPUT 

C  NUM=SAMPLE  SIZE 
C  FAIL=VECT0R  OF  SAMPLE  DATA 
C  OUTPUTS 

C  ESTIMATES  OF  LOCATION  PARAMETER  CW«  SHAPE  PARAMETER 
C  KU •  AM D  SCALE  PARAMETER  THETAU 

c  variables: 

C  INTEGER 
C  NUM=SAMPLE  SI7E 
C  REAL 

C  FAIL  =VECTOR  OF  SAMPLE  DATA 
C  RANK=VECTOR  OF  MEDIAN  RANKS 
C  F  A I L  MC  =  FAILURE  TIMES  MINUS  C 
C  LNFA I_=VECTOR  OF  LN  OF  SAMPLE  DATA 
C  CWtKWt  THET AU=  AS  ABOVE 

C  SUMl=\iUM*SUM  OF  FAILURE  TIMES**K  *LN<FAILURE  TIMES) 

C  SUM2  =S UM  OF  FAILURE  TIMES**K 
C  TOTALX=SUM  OF  LNCFAILURE  TIMES) 

C  SUMO  =S  UM1/SUM2  -TOTALX 

C**** ************************ *********** 

C 

INTEGER  NUM 

REAL  FAILMC(100)*SUM1  *SUM2*  SUMO 
REAL  FAIL(IOO)  »R  ANK  (  2  ) 

REAL  LNFAIL(10G)*CW «KU*EK*THETAU« TO TALX *X 
C  DECLARATIONS  COMPLETE 

C**************************** 

c 

D">  10  1=1*2 

RAMK<I)=C I-.3)/CNUM*.4 ) 

10  CONTINUE 

C  ESTIMATE  C  BY  LINEAR  EXTRAPOLATION  FROM  THE  FIRST 
C  TWO  ORDER  STATISTICS 

CM=FAILC2)-(FAIL(2)-FAIL(1) ) *R ANK < 2 ) /< RANK < 2>-R AN* < 1 ) ) 
C  CALCU-ATE  THE  LOG  OF  THE  FAILURE  TIMES  MINUS  C 
C  AND  HE  TOTAL  LOG  OF  THE  FAILURES  MINUS  C 
T0TALX=0. 

DO  15  I  =1  *NUM 

FAILHCII)=FAILCI)“CU 

LNF AIL( I)=ALOG(FAILMC(I )) 

TOT ALX=TOTALX*LNFAILCI) 

15  CONTINUE 
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C  CALCULATE  K  BY  AM  ITERATIVE  TECHNIQUE  OF  MAXIMUM 
C  LIKELIHOOD. 

KV=1. 

EK=2. 

90  IFC A8S(KV-EK>.LT.. 00001)60  TO  80 
SUMI=0. 

SUN2=0. 

00  20  I=1*NUM 
X=F AILMCC I)**Ky 
S  UM I  =  SUM1+LNF AIL  C I ) *X 
!  SUM2  =  SUM2*X 

|  20  COM  TINUE 

SUH0=< t SUM1*REAL(NUM> ) /SUM2 l-TOTALX 
EK=  REAL (NUM)/SUM0 
KW=  <2.*EK«-KU)/3. 

GO  TO  90 

C  CALCULATE  THETA  USING  THE  VALUE  OF  K 

C  AND  THE  SUMS  DETERMINED 

80  THETAW=CSUM2/REALCNUM)1**<1./Ky» 

RETURN 
E  <0 


PROGRAM  PAR  AM 
C  INPUT 

C  N=SAM=»LE  SIZE  (BEFORE  CENSORING)*  N=100  OR  LESS  AS 
C  DIMENSIONED 

C  SS1=0  IF  SCALE  PARAMETER  THETA  IS  KNOUN 
C  SS1=1  IF  THETA  IS  TO  BE  ESTIMATED 
C  SS2=0  IF  SHAPE  PARAMETER  K  IS  KNOUN 
C  SS2= 1  IF  K  IS  TO  BE  ESTIMATED 
C  SS3=  0  IF  LOCATION  PARAMETER  C  IS  KNOUN 
C  SS3=1  IF  C  IS  TO  BE  ESTIMATED 
C  T  ( I )  =1 • TH  ORDER  STATISTIC  OF  SAMPLE  CI=1»N> 

C  M=NUM3ER  OF  OBSERVATIONS  REMAINING  AFTER  CENSORING 
C  N-M  FROM  ABOVE 

C  C(1)=INITIAL  ESTIMATE  OR  KNOUN  VALUE  OF  C 
C  THET AC 1)=INITIAL  ESTIMATE  OR  KNOUN  VALUE  OF  THETA 
C  EK<1 )= INIT IAL  ESTIMATE  OR  KNOUN  VALUE  OF  K 
C  MR=NUMBER  OF  OBSERVATIONS  CENSORED  FROM  BELOU* 

C  NORMALLY  u  INITIALLY 
C  ********  ******* *******  ********** 

C  OUT9UT 

C  N.SS1.  SS2*SS3* M*CC1)* THETA (l).EK(l). MR 
C  --SAME  AS  FOR  INPUT 

C  CCd) =ESTIMATE  AFTER  J-l  ITERATIONS 
C  (OR  KNOUN  VALUE  OF  C) 

C  THET A(  J»  =EST IMATE  AFTER  J-l  ITERATIONS 
C  (OR  KNOUN  VALUE  OF  THETA) 

C  EK ( J )  =  E STI MATE  AFTER  J-l  ITERATIONS 
C  (OR  KNOUN  VALUE  OF  K) 

C  MAXIMUM  VALUE  OF  J  DIMENSIONED  IS  550 

C  EL=NAT  UR  AL  LOG  OF  LIKELIHOOD  FOR  C(J)*  THETA(J)*  EK(J) 
C  REFERENCE 

C  HARTER*  H. LEON  AND  A. H. MOORE.  MAXIMUM  LIKELIHOOD 
C  ESTIMATORS  OF  THE  PARAMETERS  OF  GAMMA  ANO  UEI8ULL 
C  POPULATIONS  FROM  CGMPLETE  AND  FROM  CENSOREO  SAMPLES 
C  TECHNOMETRICS,  7(1965) 


DIMENSION  T(100)*C( 550>*THET A(550 ) »EK ( 550 ) «X( 56 ) « Y ( 55) 
SSI  =1. 

SS2-1* 

SS3=1. 

N=l  0 
f»=l  0 

THET  A ( 1 ) =1 • 

C  (l  )  =  0. 

EK(  1 )  =1 . 

0 

T  (l  )  =  146.96 
T (2  >=162.52 
T  (  3  )  =  1 75.  64 
T  (4  )  =  220.447 
T (5)=223.9 
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T (6  >  =  261*63 
7(71  =  329.99 
r (8  >=334.266 
T  <9  )  =  350 . 70 
T  (l  0  >  =  359 .19 
IF(N>66, 66*104 
104  EN=N 

IF(H>66,66,110 
110  EM=  R 

31  E  LN  N=  0  • 

EMR=HR 
MRP  =  MR* 1 

33  MM=N-M*1 

90  34  I=NH,N 
EI=  I 

34  EL'JM  =  ELNM«-ALOG(EI> 

IFC  MR  >66,35*74 

74  00  75  1=1, MR 
EI=  I 

75  EL’4H=ELNM-ALOG(EI> 

35  DO  30  J=l ,550 
IF( J-l>66,25,37 

37  JJ=J-1 

j  S  K=  0  • 

:  s  u=  o . 

!  00  6  I  =  MRP,M 

|  6  SK=3K«-(T(  I  >-C(  JJ>  >**EK(  JJ> 

I F( SS 1 >7*7,8 

7  THETAC J>=THETA( JJ> 

GO  TO  9 

3  IF( MR >66,19,20 

19  THETA (J>  =  <(SK  +  (EN-EM>  *(  T(M>-C (JJ>  >**EK(JJ>  >/EM> 

1**( l./EK( JJ>> 

GO  TO  9 

20  X  <1 >=TH£T  A( JJ) 

LS=  0 

00  21  L=l*55 
LL=L-1 
LP=L*1 
X(LP>  =  X<L> 

ZRK  =  ( ( T (MRP) -C( J J>  >/X(L>>**EK(J«J) 

V(L>=-EK(JJ>* ( EM-EMR )/X(L)*EK(JJ)*SK/X(L)**(EK(JJ)»l.) 
1*E< (JJ>*(E9-EM>  * ( T ( M>-C(JJ> )••£«( JJ  >  /  X  (  L  >  *  *  (  EK  ( JJ)*1.) 
1-EMR  *EK(JJ>*ZRK*EXP(-ZRK>/(X(L>*(1.-EXP(-ZRK>>> 

IF( Y(L>  >53,73,54 

53  LS=LS-1 

I F( LS*L>38,35,38 

54  L  S=L3*1 

IF( LS-L >38,56,58 

55  X (LP>=.5*X(L> 

GO  TO  61 

56  X(LP>=1.5*X(L> 
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I 


SO  TO  61 

58  IFC  YCL)*YCLL))60«73?S9 

59  LL=  LL-1 
GO  TO  58 

60  XCL*)=XCL)»YCL)*CXCL)-XCLL) ) / CY CLL) -YC L ) ) 

61  IFC  ABSCXCLP>-XCL))-1«E-4)73»73#21 

21  CONTINUE 

73  THETACJ)=XCLP> 

9  EKC  J) =EK( JJ ) 

10  IFC  SS2)  12  tl 2*11 

11  00  17  I  =HRPf 

17  SL=SL»ALOGCTCI)-CCJJ) ) 

XCl  >  =  EKC J ) 

LS=0 

DO  51  L  =  l *55 

SL<=0. 

DO  18  I=MRP»H 

18  3LK=SLKM  ALOGCTCI  )-CC JJ ) )-ALOGC THET AC J) ) ) *C T C I ) -C C J J ) ) 
1**X  CL) 

LL=L-1 
LP=L*1 
X  CLP) =X  CL  ) 

ZRK  =  C  CT  CHRP )-CC  JJ) ) /THETACJ)  )**XCL) 

Y  CL )  =  CEM-EHR ) * C 1 • /X  CL )-ALOG C THE T A C J ) ) )*SL-SLK/ THET A C J) 
1**X  C  L  )♦•  C EN-EM)  *  C  ALOG C  THETA Cd))-ALOGCTCM)-CCJJ)))*C  TOO 
l-CC  JJ))**XCL) /THETA  CJ)  *  *XCL)*EM,R*ZRK*C  ALOGC/RK)  /XCL  )  ) 
1*EX PC -ZRK)/ Cl .-EX  PC -ZRK)) 

IFC  YCL) )43» 52*44 
A3  L5=  L3-1 

I  lFCLS+L)47»45*47 

!  44  L$=L3*1 

i  IFC L3-L)47.46,47 

45  X  CLP)=.5*XCL) 

GO  TO  50 

46  XCLP)=1.5*XCL) 

GO  TO  50 

47  IFC YCL)*YCLL) )49»52 .48 

48  L  L=  LL-1 
SO  TO  47 

49  XCLP)=XCL)*YCL) MXCD-XCLL) ) / C Y CLL) -Y CL) ) 

5  0  IFC  ABSCXCLP)-XCL))-l.E-4)52»52t51 

51  CONTINUE 

52  EKC  J) =X  CLP) 

12  CCJ)=CCJJ) 

62  IFCSS3)25«25V14 

14  IFC l.-EKC J) )16*78*78 

78  IFC  SS1*SS2)57*57«16 

16  XC1)=CCJ) 

LS=0 

DO  23  L=l»55 
SKI  =0 
3  R=  0 
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00  15  I=NRP.M 

SKI=SK1*CTCI)-XCL>)  **CEKCJ)-1.) 

15  SR=SR«-1./CTCI  )-XCL)  » 

LL=L-1 
LP=L«-1 
X  <UP)  =X  CL  ) 

ZRK  =  C CTCMRP)-XCL))/THETAC J) >**EKC J) 

V CL )=C1.-£KC J))*SR«-EKCJ)*CSK1*CEN-EM> *CTCM>-XCL>> 

1*  *CEKCJ)-1.) > /THETA (J)* *EK( JI-EMR *EKC J)*ZRK*EXPC-ZRK> 
1/ CC  T C  MRP) -X CL))*C1« -EXP  C-ZRK  )  )  > 

IFC  YCLM39.24.40 

39  LS=LS-1 

I FC  L3*L ) 7  0. 41 »70 

40  LS=LS*1 

IFC  LS-L >70.42.70 

41  X CLP)=.5*XCl> 

GO  TO  22 

42  XCLP)=.5*XCL)*.5*TC l > 

GO  T 0  22 

70  IFCYCL)*YCLL)>72»24»71 

71  LL=LL-1 
GO  TO  70 

72  XCLP)=XCL)«-YCL>  *(X(L)-X{LU)/(YaL)-T(L)» 

22  IFC  ABSCX CLP) -XCL>>-1. £-4)24,24. 23 

23  CONTINUE 

24  CCJ)=XCLP> 

GO  T0  25 

57  C(J)=T(1) 

25  IFC  MR  >66*39  *69 

38  00  63  1  =  1. M 

IFC  CC J)»l.E-4-TCI>>  68. 6  7.67 

67  MP=MR«>1 

63  CCl >=TC1> 

68  IFC *3)66.69.31 

69  SK=0. 

SL=0. 

00  36  I =MRP»M 
SK=SK+CTCI)-CCJ> >**EKCJ> 

36  SL=SL«-ALOGCTC I)-CCJ >> 

ZRK  =  C  CTCHRP>-CCJ))/IHETAC J) )**EKCJ» 

EL=ELNM-MEM-EMR)*CALOGCEKCJ))-EKC  J)  *ALOGC THETAC J)  >  >  ♦• 

1 C  E<  C  J)-l.)*SL-CSK*CEN-E>')*CTCH)-CCJn**EKCJ)>/CTHETA 
1C J) **EKCJ))*EMR*AL0GC1« -EXP  C -ZRK ) ) 

IFC  J-3) 30 *27.27 

27  IFC  ABSCCC JI-CCJJ) >-l. £-4)23,28, 30 

28  IFC ABSC THET AC JI-THETACJJ)  l-l • £-4)29.29.30 

29  IFC  ABSCEKCJ)-EKCJJ) )-l.E-4)66.6  6.30 

30  CONTINUE 

66  STOP 

EMD 
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APPENDIX  C  Derivation  of  Maximum  Likelihood  Equations 


The  two- parameter  Weibull  distribution  function  is 
f(t;e,k)  =  pr  exp  [-(|)kJ 

Assuming  a  random  sample  of  n  independent  failures  ,  T2»  •••  Tn# 
the  likelihood  function  is 

L(e,k)  =  i£1  exp 

The  natural  logarithm  of  the  likelihood  function  can  be  used  to  find 
the  maximum  likelihood  since  the  maximum  point  will  be  the  same  for  both 
the  logarithm  and  the  function.  The  random  samples  can  be  considered 
constants  for  the  likelihood  function.  If  the  partial  derivatives  of  the 
logarithm  of  the  likelihood  function  are  taken  with  respect  to  each  of 
the  two  parameters  and  set  equal  to  zero,  the  two  equations  can  be  solved 
for  the  two  unknown  parameters  by  first  solving  for  e  in  terms  of  k,  then 
substituting  and  solving  for  k. 

ln  L  =  a1  1  exp[-(ii)k]) 

1=1  e  0 

.  k 

■  i?1[ln  k  +  (k-1  )  ln  t^  -  k  In  e  -  (^)  ] 

.  .  k 

=  n(  ln  k  -  k  ln  e)  +  t  (k-1  )  ln  -  ( H)  ] 


3  ln  L  = 

3* 

0  - 


nk  k 

JT  +  ^kTT 

nk  k 

e  +  ©k+1 


n 

i£i 


set  equal  to  zero 


n 


9 


k 


lk 
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3  In  L 

nr" 


n 

k 

k 


-  2  -  n  In  •  ♦  jfj  t  In  ti  -  (^)k  ln(^)] 

=  £  -  n  In  o  +  i£1  In  (^i)k  In  (£i) 


k 

n  _  n 
k  k 


f'.z.  t  k 


In  p1  — )  +  i£iln  tA  -  n  ",  k  i|i  ln<H) 


ir=l 


n  n,/n  lAn  n 

k  •  k  ln(i?-a  )  k  ln  n  +  iE=i  ln  ti 


n 

zn  t.k 

n  .  k 

iii‘i 

1=1  i 

n  n 

r  +  -  z  , 
k  i  =  l 

111  - 

i£l  i 


+  £  ln  n  -  k . z  t . k  iti  tik  ln  n  +  k.z.t.k  i=lti  ln  i 

1  =  1  1  1=1  1 


-  £  m  .?  tik 

k  i  =  l  1 


n  n  n 
,  +  ,  In 


nth  .2, 


ti K  ln  ti 


k  i£l“  i  ’  1  =  1  l 

n 

-  ln  t±  -  n  t  Tc  ln 


i  =  l  l 
n 


n  .2  t.  ln  t. 

-izl  ■  L- - i  -  .2,  In  t. 

n  ,  k  i=l  i 


ill 


set  equal  to  zero 
and  solve 
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APPENDIX  D  Component  Failure  Data  -  Time  To  Failure 


Component  1 


139 

146 

162 

175 

220 

233 

261 

290 

308 

329 

334 

335 

350 

359 

463 

Component  2 

113 

216 

246 

252 

296 

326 

415 

431 

444 

480 

529 

535 

661 

664 

995 

Component  3 

1 63 

174 

272 

289 

305 

348 

353 

371 

384 

462 

463 

508 

510 

518 

570 

Component  4 

92 

99 

142 

150 

151 

152 

160 

167 

182 

187 

200 

201 

21 1 

214 

216 

Component  5 

160 

186 

200 

204 

206 

213 

259 

261 

265 

287 

298 

304 

326 

327 

438 

Component  6 

107 

111 

142 

162 

189 

212 

231 

233 

266 

274 

306 

369 

498 

509 

719 

79 


VITA 


Murray  Ross  MacDonald  joined  the  Royal  Canadian  Air  Force  in 
April  1964.  He  completed  navigation  training  and  received  his  wings  in 
June  1967  after  which  he  served  as  a  navigator  and  tactical  coordinator 
on  Neptune  and  Argus  long-range  patrol  aircraft.  He  completed  the 
year-long  Aerospace  Systems  Course  in  June  1973,  then  served  as  a  project 
officer  in  the  Aerospace  Engineering  Test  Establishment.  In  June  1980  he 
completed  a  Bachelor  of  Science  in  Mathematics  degree.  He  was  employed 
in  NORAD  Headquarters  as  the  Branch  Head  for  deep  space,  then  the 
Division  Chief  for  the  Space  Detection  and  Tracking  network  until 
entering  the  School  of  Engineering,  Air  Force  Institute  of  Technology,  in 
June  1 981 . 


80 


SECURITY  q  ^SSlHi.ATlQN  OF  THIS  PAOfc  (When  liutefintated) 


13  number  OF  PACES 
88 _ 


IS*.  DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 

1ft.  DISTRIBUTION  STATEMENT  (ot  it it*  KeporfJ 

Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (ot  the  mbetrect  entered  In  Block  20.  II  different  from i  Report) 


19-  KEY  WORDS  f  Confirm*  on  tmvoreo  elde  It  necoeemry  and  Identify  by  block  number) 

Monte  Carlo 

Reliability  Confidence  Limits 

Weibull  Distribution 

20.  ABSTRACT  ( Continue  on  reveree  »ldo  It  n«c«n«ry  and  Identity  by  block  number) 

A  Double  Monte  Carlo  method  of  obtaining  confidence  limits  for  complex 
systems  based  on  component  failure  data  assuming  a  three  parameter  Weibull 
distribution  was  developed.  Three  new  parameter  estimation  routines  were 
developed  and  compared  with  the  Harter-Moore  three  parameter  maximum  likelihood 
routine  for  use  with  the  Monte  Carlo  method.  The  sensitivity  of  the  method  to 
system  reliability,  sample  size,  and  number  of  points  in  the  component 
reliability  distributions  was  assessed.  An  approximate  method  of  calculating 
and  correcting  for  parameter  estimation  bias  was  developed  and  illustrated. 

DO  1  JARM7I  1473  COITION  OF  t  NOV  6ft  IS  OBSOL  ETE 


15-  SECURITY  CLASS,  (o I  thta  report) 

Unclassified 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


J  RECIPIENT'S  CAT  ALOG  NUMBER 


5  TYPE  OP  REPORT  ft  PERIOD  COVERED 


MS  Thesis 


6  PERFORMING  O^G.  REPORT  NUMBER 


•  CONTRACT  OR  GRANT  NUMBERS 


10  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  ft  WORK  UNIT  NUMBERS 


12.  REPORT  DATE 

December  1982 


REPORT  DOCUMENTATION  PAGE 


REPORT  NUMBER 

AFIT/GSO, MA  82D-1 


ft  title  (end  Subtitle) 

A  MONTE  OAHU)  TECHNIQUE  SUITABLE  EUK  OBTAINING 
COMPLEX  SPACE  SYSTEM  RELIABILITY  CONFIDENCE 
LIMITS  FKOM  COMPONENT  TEST  DATA  WITH  THHEE 
UNKNOWN  P 


7  Au  T HORf  *) 


Murray  H.  MacDonald 


9  performing  organization  name  and  adoress 

Alt  Force  Institute  of  Technology  ( AFIT-EN ) 
Wr ight-Hatterson  AFH ,  Ohio  45433 


n  controlling  office  name  ano  aooress 


ft  MONITORING  AGENCY  NAME  ft  ADDRESSfif  different  from  Controlling  Office) 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  fWhtfl  Doto  Entered) 


^tCUW.TY  CLASSIFICATION  OF  Tn.S  P»QE(Wh,,  p...  Enl.cd. 


from  74%  to  96%  wi  th1Lmp<t,nentdfaUuretsa^e1effeCtive  at  system  reiiabilities 
Linear  Least  Squares  parameter  estimation  Jout^e'^elopeJ!  fiVe  *ith  the 


SECu*)Ty  CLAUiriCATlOH  or  ru**  |»A6£rirh#f»  £nUf«rf; 


